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SECTION  I 


INTRODUCTION 

The  purpose  of  this  study  is  to  certify  the  ability  of  a  particular 
timizatior  code  to  solve  trajectory  optimization  problems  associated  with 
hypervelocity  vehicles,  that  is,  the  maximum  crossrange  problem  and  the 
maximum  plane  change  problem.  The  optimization  problems  are  converted  into 
parameter  optimization  problems  (nonlinear  programming  problems)  and  are 
solved  by  the  augmented-Lagrangian  method.  The  particular  optimization  code 
which  will  be  used  has  been  created  at  the  Atomic  Energy  Research 
Establishment  in  Harwell,  England. 

The  physical  model  used  for  the  trajectory  problems  is  described  in 
Section  II.  In  Section  III,  the  optimal  control  problem  is  converted  into 
a  parameter  optimization  problem,  and  the  optimization  code  is  discussed 
in  Section  IV.  That  part  of  the  code  to  be  provided  by  the  user  is 
presented  in  Sec' ion  V.  The  solution  of  the  optimization  problem  is 
carried  out  in  Section  VI,  and  the  conclusions  are  presented  in 
Section  VII.  Finally,  the  documentation  provided  with  the  optimization 
code,  the  listing  of  the  optimization  code  and  user  code,  and  the 
definition  of  the  variables  in  the  user  code  is  presented  in  the 
appendices. 
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SECTION  II 


PHYSICAL  MODEL 

In  this  section,  the  physical  model  used  for  the  reentry  and  plane  change 
problems  is  discussed.  The  model  includes  the  equations  of  motion,  the  earth, 
the  atmosphere,  the  aerodynamics,  the  propulsion,  the  performance  Indices,  the 
boundary  conditions,  and  the  Inequality  constraints. 

2.1  Equations  of  Motion 

The  equations  of  motion  used  for  the  reentry  and  plane  change  problems  are 
those  for  thrusting  flight  over  a  rotating  earth.  If  the  vehicle  is  assumed  to 
be  flying  west  to  east  and  if  a  positive  bank  angle  is  assumed  to  generate  a 
heading  toward  the  north,  these  equations  are  given  by 

0  •  V  cos  y  cos  i/i/r  cos  $ 

$  ■  V  cos  y  sin  *fr/r 

r  ■  V  sin  y  (1) 

V  -  (T  cos  a  -  D)/m  -  g  sin  y  +  u>zr  cos  $  (sin  y  cos  <j>  -  cos  y  sin  $  sin  ip) 

Y  ■  (T  sin  a  +  L)  cos  y/mV  +  (V2/r  -  g)  cos  y/V  +  2u  cos  $  cos  ip 

+  (w2r/V)  cos  4  (cos  y  cos  $  +  sin  y  sin  $  cos  ip) 

p  ■  (T  sin  a  +  L)  sin  y/mV  cos  y  -  (V/r)  cos  y  cob  ip  tan  p 

+  2u  (tan  y  cos  p  sin  ip  -  sin  p)  -  (u>2r/V  cos  y)  »in  $  cos  p  cos  \p 

In  these  equations,  e  is  the  longitude,  4>  is  the  latitude,  r  is  the  radial 

distance  from  the  center  of  the  earth  to  the  vehicle  center  of  gravity,  V  is 
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Che  velocity  relative  to  the  earth,  y  is  the  flight  path  inclination,  ip 
is  the  heading  angle,  m  is  the  mass,  T  is  the  thrust,  D  is  the  drag,  L  is 
the  lift,  a  is  the  angle  of  attack,  u  is  the  bank  angle,  and  u)  is  the 
angular  velocity  of  the  earth. 

2.2  Earth 

The  earth  is  assumed  to  be  a  sphere  whose  radius  represents  mean  sea  level 
and  is  denoted  by  rg  .  As  a  consequence,  the  acceleration  of  gravity  is  given 
by  the  inverse-square  law 


g  -  gg(rs/r)2  (2) 

where  gg  is  the  acceleration  of  gravity  at  sea  level.  Also  the  altitude  of 
the  vehicle  above  mean  sea  level  satisfies  the  relation 

h  -  r  -  rg  (3) 

The  values  of  the  constants  associated  with  the  characteristics  of  the  earth 
are  listed  below: 

rg  -  20926428  ft 

gg  -  32.174  ft/sec2  (4) 

u  ■  7.2921151E-5  rad/sec 

2. 3  Atmosphere 

In  order  to  obtain  atmospheric  properties  as  a  function  of  altitude,  the 
approximate  atmosphere  known  as  the  1962  U.S.  Standard  Atmosphere  has  been 
employed,  with  the  additional  assumption  that  the  composition  of  the  atmosphere 
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is  constant.  Here,  the  atmosphere  is  divided  into  a  number  of  layers  in 
which  the  temperature  gradient  is  assumed  constant.  Bence,  the  absolute 
temperature  in  each  layer  is  given  by 

t  -  t  +  1  (h  -  h  )  ,  h  >  h  (5) 

n  n  n  —  n 

where  is  the  absolute  temperature  at  the  beginning  of  the  layer,  1r 

is  the  constant  temperature  gradient,  and  h  is  the  altitude  at  the  beginning 

n 

of  the  layer.  Once  the  temperature  is  known,  the  speed  of  sound  is  obtained 
from  the  relation 

a  -  <kR  t)1/2  (6) 


where  k  is  the  ratio  of  specific  heats  of  air  and  R  is  the  gas  constant  of 
air  at  sea  level.  For  those  layers  where  the  temperature  gradient  is  nonzero, 
the  density  ratio,  o  -  p/pg  ,  can  be  expressed  as 


o 


C 

n 


-(1  +  gQ/R  1  ) 
t  an 


(7) 


where  Pc  is  the  density  at  sea  level  and  C  is  a  known  constant  for  each 
°  n 

layer.  For  those  layers  where  the  temperature  gradient  is  zero,  the  density 
ratio  becomes 


o  ■  C  exp  (-g  h/R  t  )  (8) 

n  on 

The  values  of  the  constants  associated  with  the  atmosphere  are  presented 
below  and  in  Table  1: 

k  -  1.4 

R  -  3086.9629  ft2/sec2  *K 

Ps  -  2.3769E-3  slugs/ft3  (9) 
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Table  1.  Atmospheric  Constants 


Layer 

(n) 

hn 

(ft) 

Tn 

(°K) 

^■n 

(°K/ft) 

C 

n 

1 

0 

288.15 

-1.9812E-3 

3.401824655257E-11 

2 

36,089 

216.65 

0 

1. 683376997149E+00 

3 

65,617 

216.65 

3.0480E-4 

9. 817858914969E+80 

4 

104,987 

228.65 

8. 5344E-4 

1 . 506414967722E+29 

5 

154,199 

270.65 

0 

4 . 394749884481E-91 

6 

170,604 

270.55 

-6. 0960E-4 

4.717851690435E-43 

”7 

/ 

200,131 

252.65 

-1.2192E-3 

1 . 562793740651E-22 

8 

259,186 

180.65 

0 

5 . 028946984109E+01 

During  the  optimization  process,  trajectories  can  be  obtained  which  go  to 
very  high  or  very  low  altitudes.  To  avoid  a  fatal  error  associated  with  expo¬ 
nential  overflow,  the  following  additional  features  have  been  incorporated 
in  the  model  atmosphere: 


h  1  0  :  t  -  288.15  -  1.9812E-3  h 

a  -  1.  -  2.9262913E-5  h 

h  >_  750,000:  t  -  180.65 

o  -  8.111816E-18 


(10) 


Hence,  below  sea  level,  the  density  ratio  is  assumed  to  vary  linearly  with 
altitude,  and  at  high  altitudes,  it  is  assumed  constant  at  the  value  for 


750,000  ft. 
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2.4  Aerodynamics 


The  drag  and  the  lift  are  related  to  the  drag  coefficient  Cp  and  the 
lift  coefficient  C  as  follows: 

L 


D  -  (1/2)  0SV2Cr 


L  -  (1/2)  PSV2CL 

where  S  Is  the  aerodynamic  reference  area.  The  drag  and  lift  coefficients 
can  be  expressed  in  terms  of  the  axial  and  normal  force  of  coefficients  as 


C_  *  C.  cos  a  +  C.,  sin  a 
DA  N 


*  CN  cos  a  -  CA  sin  o 

The  axial  force  coefficient  is  composed  of  a  skin-friction  term  and  a  pressure 
term,  that  is, 


C.  -  C.  +  C. 
A  ASF 


where 


C,  -  A.h2  +  B.h  +  C, 
ASF  1  11 


CA_„  “  V  +  V  +  C2 


In  these  relations,  A^  ,  B^  ,  and  are  known  functions  of  the  Mach 

number,  M  ■  V/a  ,  while  k^  ,  and  are  known  functions  of  Mach 

number  and  angle  of  attack.  Finally,  the  normal  force  coefficient  is  given 
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CN  "  A3“2  +  B3°  +  C3 


(15) 


where  A^  ,  ,  and  are  known  functions  of  Mach  number. 

The  values  of  the  constants  associated  with  the  aerodynamics  are  given 
below  and  in  Tables  2  through  4. 


S  -  125.84  ft2  (16) 


Table  2.  Axial  Skin-Friction  Force  Constants 


M 

A1 

B1 

C1 

0.2 

0 

7.80R-3 

.0114 

1.2 

0 

7.70E-8 

.0076 

5.0 

0 

5.60E-8 

.0028 

10. 

2.40E-12 

-7.11E-7 

.0578 

i  20- 

1.38E-12 

-3.81E-7 

.0297 

Table  3.  Axial 

Pressure  Force  Constants 

M 

A2 

B2 

C2 

> 

a  <  16® 

a  >  16° 

a  <  16° 

a  >  16® 

a  <  16° 

“a  >.  16 4 

0.2 

-1.00E-4 

-1.00E-4 

-2.19E-3 

-2.19E-3 

.0350 

.0350 

1.2 

-7.81E-5 

-7.81E-5 

-1.25E-3 

-1.25E-3 

.0760 

.0760 

5.0 

1.09E-5 

4.86E-6 

-9.62E-4 

-1.13E-4 

.0324 

.0204 

10. 

1.41K-5 

-6.60E-6 

-9.00E-4 

3.49E-4 

.0261 

.0114 

>  20. 

1.33E-5 

-6.25E-6 

-8.19E-4 

3.58E-4 

.0243 

.0105 
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Table  4.  Normal  Force  Constants 


M 

A3 

B3 

C3 

0.2 

1.33E-4 

2.68E-2 

-.030 

1.2 

-7.81E'-5 

2.90E-2 

-.030 

5.0 

2.94E-4 

1.41E-2 

-.035 

10. 

4.38E-4 

6.50E-2 

-.035 

>  20. 

4.38E-4 

6.50E-2 

-.035 

In  the  computer  program,  the  above  tables  are  read  by  linear  Interpolation. 
Also,  for  N  >  20  ,  the  M  ■  20  value  Is  used  because  the  extrapolation  of  the 
above  numbers  could  lead  to  wrong  results. 

2.5  Propulsion 

For  the  plane  change  problem,  rocket  engines  are  used  to  aid  In  the  plane 
change  and  to  increase  the  energy  of  the  vehicle.  It  is  assumed  that  the  engines 
bum  once  and  that  the  propellant  mass  flow  rate  and  the  thrust  are  constant 
during  the  bum.  Hence,  during  the  bum,  the  mass  of  the  vehicle  satisfies 
the  relation 


m  ■  m^  -  8(t  -  tfc) 


(17) 


where  B  is  the  mass  flow  rate  and  where  and  t^  are  the  mass  and  the 

time  at  the  Initiation  of  the  burn. 

The  values  of  the  constants  associated  with  the  propulsion  characteristics 
are  as  follows: 
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m^  ■  nig  ■  335.67477  slugs 
6  ■  0.34768573  slugs/sec 

T  -  3300  lb  (18) 


m  ■  179.18195  slugs 
P 

At,  ■  515.35606  sec 

D 

where  m  is  the  propellant  mass  and  At.  is  the  total  bum  time.  In  the 
p  b 

reentry  program,  the  propellant  mass  is  assumed  to  have  been  expended,  so  the 
initial  mass  is  m^  *  156.49282  slugs. 

2.6  Performance  Indices 

In  the  reentry  problem,  it  is  desired  to  maximize  the  crossrange.  Hence, 
to  have  a  minimization  problem,  the  performance  index  is  written  as 


J  -  -  *f  (19) 

For  the  orbital  plane  change  problem,  it  is  desired  to  maximize  the 
plane  change.  Maximizing  the  plane  change  is  equivalent  to  maximizing  the 
orbit  Inclination  because  the  initial  orbit  inclination  is  zero.  The  orbit 
inclination  is  the  angle  between  the  Inertial  angular  momentum  vector 


H 


T"  X  V 

r  inertial 


(20) 


and  the  angular  velocity  vector  of  the  earth,  where  the  bar  denotes  a 
vector.  From  the  definition  of  the  dot  product,  it  is  seen  that 


cos  i  ■  (H  *  ii))/Hu 


(21) 
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Also,  the  inertial  velocity  is  given  by 


Vinertial 


V  +  w  x  r 


(22) 


Finally,  If  the  vector  operations  are  carried  out,  the  problem  of  maximizing  i 
becomes  that  of  minimizing  the  cosine  of  the  Inclination  at  the  final  point: 


J  ■  (cos  i)^ 

where 


(23) 


cos  1 


cos  $  (V  cos  Y  cos  +  rv>  cos  $) 


[(V  cos  y  cos  <J>  +  ra>  cos  $)2  +  (V  cos  V  sin  >J02] 


T 71 


(24) 


2.7  Prescribed  Boundary  Conditions 

The  Initial  conditions  for  both  problems  are  as  follows: 

tQ  -  0  ,  eQ  -  0  ,  -  0  ,  hQ  -  364,800 

VQ  -  24,500  ft/sec  ,  yQ  -  -1.2  deg  ,  g/Q  «  0 

For  the  reentry  problem,  the  final  condition  is  given  by 

hf  -  100,000  ft 

The  final  conditions  for  the  plane  change  problem  are  the  following: 
h^  ■  608,000  ft  ,  J  23,500  ft/sec  ,  -  0 


(25) 


(26) 


(27) 
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2.8  Inequality  Constraints 


The  angle  of  attack  and  the  load  factor ,  n  -  L/W  ,  are  required  to 
satisfy  the  Inequality  constraints 

a  ^  40  deg  ,  n  ■  _<  4,5  (28) 


2.9  Nondimens ional  Variables 

Nondimensionalizatlon  has  the  effect  of  scaling  the  variables  and  can  be 
very  beneficial  to  optimization.  In  this  system,  all  angles  are  in  radians. 
Furthermore,  to  simplify  the  notation,  a  nondimens ional  variable  is  denoted 
by  the  same  symbol  as  the  dimensional  variable  and  a  tilde.  Hence,  the 
following  nondimens ional  variables  are  defined: 

(29) 
=  B 


t/(rs/8s)1/2  "  *  *  v/<Vs)1/2  =  ^  *  r/rs  s  *  »  ajtrs/8s)1/2  E  “ 

2/2 

T/mo8s  5  T  ’  D/no8s  -  D  •  L/mogs  5  L  ’  “/no  5  m  *  8(rs/8s)  /mo 


M 
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SECTION  III 


FORMULATION  OF  THP  OPTIMIZATION  PROBLEM 

In  this  section,  the  optimal  control  problem  is  stated  in  terms  of  the 
nondimensional  variables,  and  the  parameter  optimization  problem  is  formulated. 

3.1  Optimal  Control  Problem 

The  optimal  control  problem  is  to  find  the  angle  of  attack  and  roll 
angle  histories  which  maximize  the  crossrange  for  the  reentry  problem 
which  maximize  the  orbit  inclination  (or  minimize  its  cosine)  for  the  plane 
change  problem.  In  terms  of  the  nondimensional  variables  presented  in 
Equation  (29),  the  optimization  problem  can  be  stated  as  the  following 
minimization  problem: 

Find  the  control  variable  histories  a(t)  and  \j(t)  which  minimize 
the  performance  index 


subject  to  the  differential  constraints 

de/dt  ■  V  cos  y  cos  5>/ ( r  cos  $) 
d$/dt  -  V  cob  y  sin  Ji/r 
dr/dt  “  V  sin  y 

dV/dl  ”  (T  cos  a  -  6)/i- sin  y/rz  +  u2x  cob  |(sln  y  cos  $ 
-  cos  y  sin  I  Bin  iji) 
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.  ..MX-.. 


dy/dl  *  (T  sin  A  +  L)  cos  u/mV  +  V  cos  y/r  -  co?  y/r2V  +  2w  cos  <j>  cos  i 
+  (w^r/V)  cos  $(cos  y  cos  $  +  sin  y  sin  $  sin 

(31) 

dji/dt  ■  (T  sin  a  +  L)  sin  y/ mV  cos  y  -  (V  cos  y/r)  cos  ij>  tan  $ 

+  2w(tan  y  cos  $  sin  -  sin  $) 

-  (ai2r/\7  COS  y)  sin  |  cos  £  cos  ^  , 


subject  to  the  prescribed  boundary  conditions 


tn  =  0  ,  Gn  -  0  ,  «L  -  0  ,  r  »  1.017432502 


V  »  0.94420438  .  yft  >■  -2 .0943951E-2  ,  k  -  0 


Reentry: 


r  =  1.0047786464 


Plane  Change:  rf  =  1.029054170  ,  Vf  *  0.90566542  ,  yf  =  0 


and  subject  to  the  inequality  constraints 


a  <_  0.69813170  ,  n  <.  4.5. 


For  the  reentry  problem,  the  vehicle  is  gliding  with  the  engine  off.  Hence, 
the  thrust,  the  mass  flow  rate,  and  the  mass  satisfy  the  relations 


T  ■  0  ,  8  ■  0  ,  m-1 


while  the  drag  and  lift  are  given  by 


D  -  15268.635  C^oV^ 


L  -  15268.635  C  oV2 

L 
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For  the  plane  change  problem,  the  engines  are  ignited  at  ,  whose  optimal 
value  is  to  be  determined,  and  burn  for  the  time  period  At,  ■  0.63901697  . 

D 

Here,  the  thrust,  the  mass  flow  rate,  and  the  mass  are  given  by 


e  - 


0  engine  off 

10.83533982  engine  cn 


m 


l.  -  S(t  -  ?b>  , 

0.46620367 


0  <  t  <  tu 

D 

So  -S  ^  1 


So 


(36) 


while  the  drag  and  the  lift  become 


D  -  932.34367  C^aV2 
l  -  932.34367  C^2 


(37) 


The  reason  for  the  different  constants  in  Equations  (36)  and  (37)  is  that  the 
initial  mass  for  the  reentry  problem  is  less  than  the  initial  mass  for  the 
plane  change  problem  (see  Section  2.5). 

3 . 2  Parameter  Optimization  Problem 

In  general,  the  optimal  control  problem  of  the  previous  section  can  be 
stated  in  standard  matrix  notation  (Reference  1)  as  follows: 

Mlmimize  the  performance  index 

J  •  <t(tf,  xf)  (38) 
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subject  to  the  differential  constraints 


x  -  f(t,  x,  u)  , 

the  prescribed  boundary  conditions 

tQ  -  0  ,  xQ  =  given  ,  y(tf,  xf)  -  0  , 
the  control  variable  Inequality  constraints 


(39) 


(40) 


C(t,  x,  u)  _>  0  ,  \  (41) 

and  the  state  variable  Inequality  constraints 

S(t,  x)  >  0  .  (42) 

In  the  parameter  optimization  method  which  will  be  used  to  solve  these 
problems,  the  inequality  co  .  • aints  must  be  in  an  algebraic  form.  This  can 
be  accomplished  by  forming  an  integral  over  the  region  where  the  inequality 
constraint  is  violated.  In  other  words,  if  e(t)  >_  0  is  a  scalar  inequality 
constraint,  an  appropriate  integral  constraint  is 

,cf  - 

E  ■  -  /  [min  (e,  0)]2  dt  >_  0  (43> 

C0 

This  integral  constraint  can  be  converted  into  a  differential  equation  and  an 
algebraic  inequality  constraint.  Let 

*n+l  ”  *(min  (e’  0))2  •  xn+l , 0  ”  0  (44) 

so  that  Equation  (43)  becomes 


E 


x  .  ,  >  0 

n+1 , f  — 


(45) 
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If  the  inequality  constraints  are  converted  into  differential  equations 
and  algebraic  constraints,  the  optimal  control  problems  can  be  restated  as 
follows: 


Minimise  the  performance  index 

J  ■  * (tf,  yf) 

subject  to 

y  -  f(t,  y,  u) 

C0  “  0  ’  y0  '  glven» 


(46) 


(47) 


¥  (tf.  yf)  -  o  ,  (48) 

and  to 

Q(tf,  yf)  >  0  ,  (49) 

The  vector  y  now  contains  the  original  state  x  and  the  states  introduced 
by  converting  the  inequality  constraints. 

The  conversion  of  this  optimal  control  problem  into  a  parameter  optimization 
or  nonlinear  programming  problem  is  performed  in  two  steps.  The  first  step  is 
to  convert  the  final  time  into  a  parameter,  and  the  second  step  ic  to  para¬ 
meterize  the  control  histories. 

The  final  time  can  be  made  into  a  parameter  by  introducing  the  transformation 

t  -  t/tf  .  (50) 
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This  transformation  allows  the  optimal  control  problem  to  be  rewritten  as 
Minimize  the  performance  index 

J  -♦  (yf.  tf)  (51) 

subject  to 


to 


y*  •  g(t.  y.  u,  tf) 

t0-0  .  yQ  =  given  ,  Tf  -  1  , 


(52) 


f  (yf,  t£)  -  0  ,  (53) 

and  to 

0(yf.  t{)  >  0  .  (54) 

Note  that  the  prime  denotes  a  derivative  with  respect  to  t  and  that  this 
transformation  fixes  the  interval  of  integration. 

The  control  histories  can  be  parameterized  in  a  number  of  ways,  but  the 
simplest  is  to  use  a  set  of  nodal  points  and  create  the  function  by  linear 
interpolation  (see  Figure  1  where  a  and  u  are  control  variables).  This 
approach  allows  the  ktJl  control  to  be  written  in  the  form 

\  -  UjJt,  ak)  (55) 

where  a^  is  a  vector  of  parameters,  that  is,  the  nodal  points  defining  the 
control  history  . 

At  this  point,  the  unknown  parameters  are  lumped  into  a  single  vector 
X  defined  as 
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X  ■  [t^,  •  •  •  »  (56) 

where  p  is  the  number  of  control  functions.  For  the  plane  change  problem, 

X  also  contains  the  Ignition  time,  ,  as  the  last  parameter.  The  total 
number  of  parameters  Is  N  .  If  values  are  given  to  the  components  of  X  , 
that  Is,  If  the  final  time,  the  control  histories,  and  the  ignition  time  are 
known,  the  differential  equations  (52)  cen  be  integrated  from  tq  ■  0  to 
-  1  to  form  the  function  ■  yf(X)  .  In  turn,  can  be  eliminated 

from  the  performance  Index  and  the  constraints  to  form  the  parameter  version 
of  the  optimal  control  problem: 

Minimize  the  performance  index 

J  -  F(X)  (57) 

subject  to  the  equality  constraints 

C±(X)  -  0  ,  1  -  1  K  (58) 

and  the  inequality  constraints 

Ct(X)  >0,  1  •  K  +  1  -»  M  (59) 
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SECTION  IV 


PARAMETER  OPTIMIZATION  METHOD 


In  this  section,  the  theoretical  basis  of  the  method  which  has  been 
■elected  for  solving  the  optimization  problem  is  discussed.  Also,  some 
remarks  are  made  about  the  code. 


4.1  Augmented- Lagranglan  Method 

The  augmented-Lagrangian  method  is  a  particular  formulation  of  the  penalty 
function  method  which  allows  convergence  to  be  achieved  without  having  to  drive 
the  weights  to  infinity  (see  Reference  2).  Here,  the  performance  index  is 
defined  as 

m 

j  -  pfx)  +  (i/2)  l  o.  r<2(x,  e.)  (60) 

i-l  1  1  1 


where 


CA(X)  -  flt 


,  i  -  1  ■*  K 


r±&,  e1) 


(61) 


min  [Ci(X)  -  9±t  01  ,  i  ^  K  +  1  ♦  H 


and  where  is  a  positive,  constant  weight. 

A  description  of  the  computational  algorithm  is  as  follows: 


1.  Guess  initial  values  for  X  ,  ,  and  a ^  -  Since  x  represents 

the  final  time,  the  controls,  and  perhaps  the  ignition  time,  it  should 
be  easy  to  run  some  constant  control  trajectories  and  find  a  value  for 
X  which  gives  a  low  value  to  the  performance  index  P  .  For  the 
first  run,  the  constants  6^  are  usually  oet  equal  to  zero.  Finally, 
the  weights  can  be  obtained  by  requiring  each  constraint  term  in 

(60)  to  have  the  same  order  of  magnitude  as  F  . 
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2.  Minimize  the  performance  index  with  respect  to  X  holding  and 

constant.  This  step  represents  an  unconstrained  minimization 
which  is  performed  wl th  the  Davidon-Fletcher-Powell  method. 

3.  Assuming  that  X  -  X(0)  ,  minimize  the  performance  index  with  res- 

I 

pect  to  0  ,  starting  from  the  X  which  results  from  step  (2). 

The  purpose  of  this  step  is  to  find  new  values  of  6  such  that  the 
constraints  are  satisfied.  Only  a  single  iteration  of  the  optimization 
is  performed.  For  equality  constraints,  this  is  done  with  the 
Newton-Raphson  method,  while  for  inequality  constraints,  a  minimi¬ 
zation  problem  is  formed  which  guarantees  that  the  corresponding 
Lagrange  multipliers  remain  nonnegative. 

4.  If  the  rate  at  which  a  constraint  is  being  reduced  is  not  sufficiently 
large,  the  value  of  the  corresponding  weight  0^  is  Increased. 

5.  If  convergence  has  not  been  achieved,  go  to  2. 

4 . 2  Optimization  Cone 

The  code  which  '^s  been  selected  for  solving  the  optimization  problem  has 
been  obtained  from  the  AFRE  in  Harwell,  England.  It  is  a  general  purpose  code 
for  solving  the  nonlinear  programming  problem  using  the  augmented-Lagranglan 
method.  The  code  contains  five  subroutines  which  are  described  briefly  below. 

VF01A:  This  subroutine  directs  the  optimization  process. 

VA09A:  This  subroutine  performs  unconstrained  minimization  using  the 

Davidon-Fletcher-Powell  method.  The  metric  is  represented  in  a 
a  factored  form  to  preserve  positive  definiteness  so  that  a  crude 
one-dimensional  search  can  be  used.  Convergence  is  superlinear. 
MC11A:  This  set  of  subroutines  (MC11A,  B,  C,  D,  and  E)  performs  matrix 
manipulations  such  as  factorization  and  the  metric  update. 
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VE04A:  This  subroutine  performs  the  optimization  for  satisfying  the 
constraints. 

VF01Z:  This  subroutine  ferns  the  augmented  Lagrangian  (60)  and  its 

derivative  from  the  performance  index  (57) ,  the  equality  constraints 
(58),  the  inequality  constraints  (59),  and  their  derivatives. 

The  documentation  which  has  been  supplied  with  these  subroutines  is  contained 
in  Appendix  A.  Also,  a  listing  of  the  code  is  presented  in  Appendix  B.  That 
part  of  the  code  which  must  be  supplied  by  the  user  is  shovn  in  Appendix  C. 

The  user  must  provide  a  subroutine  entitled  VF01B  which  computes  the 
performance  index  F  and  the  constraints  C  as  well  as  the  derivatives  F^ 
and  .  In  order  to  minimize  the  amount  of  set-up  time,  it  has  been 
decided  to  compute  the  derivatives  numerically,  using  central  differences. 

This  approach  also  allows  the  user  to  change  the  physical  model  without  having 
to  change  derivative  subroutines  as  well. 

This  optimization  code  has  been  verified  by  applying  It  to  the  solution 
of  a  number  of  problems  whose  solutions  are  known.  These  problems  include  a 
simple  quadratic  function  in  two  variables  with  a  linear  equality  or  inequality 
constraint  and  the  Rosenbrock  function  with  a  quadratic  constraint.  Also,  the 
optimal  control  problem  known  as  the  lunar  launch  problem  has  been  solved. 

The  problem  Is  to  transfer  a  rocket  from  the  surface  of  the  moon  to  orbital 
conditions  in  minimum  time,  assuming  that  the  ratio  of  the  thrust  to  the  mass 
is  constant.  The  control  variable  is  the  angle  between  the  thrust  vector  and 
the  horizontal.  To  check  out  inequality  constraints,  an  active  upper  limit 
has  been  imposed  on  the  steering  angle.  In  all  cases,  the  Harwell  code  per¬ 
formed  extremely  well,  and  the  known  extremal  solutions  were  obtained  in  a 
correct  or  reasonable  number  of  iterations. 
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After  solving  the  above  problems,  it  became  apparent  that  some  minor  modi¬ 
fications  had  to  be  made  in  the  Harwell  code.  First,  if  an  inequality  constraint 
was  imposed  and  was  exactly  equal  to  zero  on  the  first  iteration,  the  program 
aborted.  To  prevent  this,  the  constraint  was  set  equal  to  a  small  positive 
value.  Second,  for  solving  more  complex  problems,  it  became  apparent  that  it 
would  be  necessary  to  interrupt  the  computation,  store  some  results  on  a  tape, 
and  restart  the  computation.  Such  modifications  have  been  made.  Finally,  for 
solving  the  reentry  and  plane  change  problems,  some  procedure  for  preventing 
the  flight  path  inclination  from  becoming  -90  deg  had  to  be  Incorporated. 

If  this  happens,  the  heading  angle  equation  blows  up,  and  the  program 
terminates.  This  situation  occurs  during  the  one-dimensional  search  associated 
with  the  unconstrained  minimization.  If  too  large  a  change  is  made  in  the 
parameters,  the  vehicle  can  end  up  in  a  vertical  dive.  The  problem  has  been 
eliminated  by  monitoring  the  flight  path  inclination  on  every  integration  step. 

If  it  gets  lower  than  -60  deg  ,  the  one-dimensional  search  is  terminated,  the 
search  stepsize  is  reduced,  and  the  search  is  restarted. 
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SECTION  V 


USER  CODE 

The  user  of  Che  optimization  code  must  provide  a  main  program  and  a 
subroutine  or  sequence  of  sut routines  in  which  the  performance  index,  the 
constraints,  and  their  derivatives  are  calculated.  These  parts  of  the  computer 
program  are  discussed  here.  Because  of  the  similarity  of  the  two  trajectory 
optimization  problems,  it  has  been  possible  to  combine  them  into  a  single 
program.  The  listing  of  the  user  code  is  presented  in  Appendix  C.  The 
variables  in  the  user  code  are  defined  in  Appendix  D. 

5.1  Main  Program 

The  main  program  MRRV  performs  two  functions.  First,  all  of  the 
parameters  needed  by  VF01A  are  defined  (see  Appendix  A).  Second,  the 
process  for  interrupting  the  computation  and  storing  intermediate  numbers 
on  tape  as  well  as  reading  intermediate  numbers  from  tape  and  restarting  the 
computation  is  controlled. 

5.2  Performance  Index.  Constraints,  and  Derivatives 

Subroutine  VF01B  calls  subroutine  SG,  which  computes  the  performance 
Index  and  the  constraints,  and  subroutine  SGX,  which  computes  the  derivatives 
of  the  performance  index  and  the  constraints. 

In  subroutine  SG,  the  tables  defining  the  control  variable  histories 
are  formed  from  the  current  values  of  the  parameters.  The  initial  conditions 
for  the  differential  equations  are  defined,  and  the  differential  equations  of 
motion  are  integrated  to  the  final  time.  From  the  final  values  of  the  states, 
the  performance  index  and  the  constraint  residuals  are  computed.  While  the 
trajectory  is  computed  In  nondimensional  variables,  the  printout  is  in  terms 
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of  dimensional  variables.  For  the  plane  change  problem,  the  engine  is  Ignited 

at  t,  and  shut  off  at  t.  +  Lt.  .  Finally,  if  the  flight  path  inclination 
b  b  d 

becomes  less  than  -60  deg  at  any  point  of  the  trajectory,  the  Integration 
is  stopped,  and  the  computation  is  redirected  to  the  one-dimensional  search 
for  a  decrease  in  the  search  stepsize. 

The  Integration  is  performed  using  a  fixed-step,  fourth-order,  Runge-Kutta 
Integration  method  (subroutine  RUNGE).  The  number  of  integration  steps  has 
been  chosen  by  making  calculations  of  the  penalized  performance  index  and  its 
derivatives  for  different  step  sizes  and  choosing  the  largest  stepsize  which 
gives  reasonable  accuracy.  This  has  been  done  to  minimize  the  computation  time 
and  to  keep  the  Integration  as  far  away  from  round-off  error  as  possible. 

The  integration  subroutine  makes  repeated  calls  to  the  subroutine  DERIV 
which  contains  the  differential  equations  of  motion  of  the  vehicle  and  the 
differential  equations  for  the  Inequality  constraints.  Here,  the  tables 
containing  the  controls  are  read;  the  aerodynamic,  propulsion,  and  mass  char¬ 
acteristics  are  determined;  and  the  right-hand  sides  of  the  differential 
equations  are  computed.  The  aerodynamic  characteristics  are  obtained  from  a 
call  to  the  subroutine  AERO  which  reads  the  aerodynamics  tables.  To  obtain 
the  atmospheric  properties, .AERO  calls  subroutine  ATM62. 

Subroutine  SLIN1  performs  linear  interpolation  of  a  table  of  data  points. 
It  has  been  provided  to  read  the  control  tables  and  the  aerodynamics  tables. 

The  last  subroutine  supplied  by  the  user  is  SGX.  Here,  the  derivatives 
of  the  performance  index  and  the  constraints  are  computed  by  central 
differences.  Hence,  two  function  evaluations  are  made  for  each  derivative 
computed. 
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SECTION  VI 


SOLUTION  OF  THE  OPTIMIZATION  PROBLEMS 

The  general  procedure  which  has  been  followed  in  the  solution  of  the 
optimization  problems  is  to  find  a  nonoptimal,  constant  control  trajectory 
which  gives  a  reasonable  value  to  the  performance  index  and  is  somewhat  near 
the  desired  final  conditions.  This  is  done  by  varying  the  parameters 
(final  time,  angle  of  attack,  bank  angle,  and  Ignition  time  if  necessary) 
systematically  until  a  reasonable  trajectory  is  obtained.  Then,  these 
parameters  are  used  as  the  initial  guess  of  the  optimal  trajectory. 

6.1  Reentry  Problem 

For  the  reentry  problem,  a  reasonable  set  of  nominal  parameters  has 
been  found  to  be  tf  -  3200  sec  ,  a  ■  20  deg  ,  and  u  ■  60  deg  .  The 
results  from  integrating  the  equations  of  motion  for  these  values  of  the 
parameters  are  presented  in  Table  5.  Note  that  the  nominal  crossrange  is 
-  26.4  deg  and  that  the  final  condition  is  given  by 

C:  -  h^/ 100,000  -  1  -  0.136  (62) 

The  main  program  contains  all  of  the  input  quantities  and  <s  presented 
in  Appendix  C  for  the  reentry  problem  (IPROB  •  1).  The  iteration  process 
starts  from  the  nominal  path  (IRS  “  0),  and  no  more  than  TMAX  *  1200  sec  of 
computer  time  are  to  be  used  on  the  first  run.  To  allow  the  program  to  reach 
the  point  where  the  time  check  is  made,  the  external  time  limit  has  been  set 
at  1500  sec.  A  total  of  eleven  parameters  (N  -  11)  is  being  used,  of  which 
one  is  the  final  time,  five  (NA  -  5)  are  the  ordinates  of  the  angle  of  attack 
table,  and  five  (NM  ■  5)  are  the  ordinates  of  the  bank  angle  table.  There  is 
one  constraint  (M  ■  1)  which  is  an  equality  constraint  (K  ■  1).  Next,  the  first 
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Table  5.  Nominal  Path  for  Maximum  Crossrange 
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parameter  X(l)  is  set  equal  to  the  nominal  final  time,  X(2)  through  X(6) 
are  set  equal  to  the  nominal  angle  of  attack  history,  and  X (7 )  through  X(ll) 
contain  the  nominal  bank  angle  history.  Also,  the  values  of  the  normalized 
time  at  each  of  the  nodes  are  the  same  for  each  table,  TTA(I)  and  TTM(I)  , 
and  are  given  by  0.0  ,  0.25  ,  0.5  ,  0.75  ,  and  1.0  .  The  unconstrained 

minimization  will  terminate  when  the  relative  change  in  each  parameter  between 
two  iterations  is  less  than  EPS  (I)  -  l.E-4  *  X(I)  .  On  the  other  hand,  the 
constraint  satisfaction  iteration  will  terminate  when  the  constraint  residual 
satisfies  the  relation  C(l)  <_  AKMIN  *  C(M  +  1)  where  AKMIN  -  l.E-4  and 
where  C(M  +  1)  is  defined  below.  Since  C(M  +  1)  “  0.136  ,  an  accuracy 
of  l.E-5  is  being  requested  in  the  constraints. 

The  expected  change  in  the  performance  index  on  the  first  iteration  is 
DFN  a  0.5  .  MAXFN  is  set  equal  to  a  large  number  (10,000)  because  the  program 
is  now  interrupted  on  a  time  basis  rather  than  a  function  evaluation  basis. 

Every  iteration  in  VF01A  is  to  be  printed  (IPR1  *  1)  as  is  every  iteration 
in  VA09A  (IPR2  “  1).  The  work  space  W(I)  is  made  to  have  IW  -  2500  words  as 
suggested  by  Appendix  A.  The  optimization  code  will  calculate  the  weight 
(MODE  »  1)  using  the  constraint  scale  factor  C(M  +  1)  *  0.136  .  To  modify 
to  code  for  interruption,  it  is  now  necessary  to  input  the  constraint  scale 
factors  in  C(M  +  I)  rather  than  in  C(I)  as  stated  in  Appendix  A.  Also, 
the  parameter  AK  -  1E60  must  be  set  here  as  well.  The  remaining  quantities 
and  the  rest  of  the  main  program  direct  the  writing  on  tape  and  reading  from 
tape  if  the  computation  is  interrupted  (CPU  time  >  TMAX) . 

A  summary  of  the  results  is  presented  in  Tables  6  though  9,  and  the 
characteristics  of  the  converged  trajectory  are  shown  in  Figures  1  through  3. 

Four  cycles  have  been  needed  to  achieve  the  desired  accuracy  in  the  final  condi¬ 
tion.  It  is  interesting  to  note  that  a  good  trajectory  is  obtained  after  one  cycle 
and  that  the  remaining  cycles  are  used  to  satisfy  the  final  condition  better.  The 
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and  h  versus 


optimal  angle  of  attack  seems  to  be  at  the  value  for  maximum  lift-to-drag 
ratio,  although  this  has  not  been  verified.  With  regard  to  the  roll  angle, 
the  vehicle  starts  out  at  a  high  roll  angle,  eo  that  it  dives  into  the  atmos¬ 
phere,  and  then  gradually  reduces  the  roll  angle  until  it  is  heading  northward 
with  a  very  small  roll  angle.  The  vehicle  achieves  a  crossrange  of  approxi¬ 
mately  -  47  deg  and  experiences  a  maximum  load  factor  of  around  1.5  g's. 

It  should  be  noted  that  the  weight  o^  has  not  been  increased  during 
the  optimization  process.  This  means  that,  after  each  cycle,  sufficient 
progress  has  been  made  in  converging  to  the  constraints.  On  the  other  hand, 
is  changed  during  each  cycle,  and  the  Lagrange  multiplier 
appears  to  be  approaching  a  limiting  value.  This  multiplier  and  the  optimal 
controls  can  be  used  to  generate  initial  values  of  the  Lagrange  multipliers 
needed  for  starting  the  shooting  method. 

The  complete  computation  required  1640  sec  of  CPU  time  (Cyber  170/750) 
and  a  total  of  3000  function  evaluations.  This  amounts  to  0.55  sec  per  function 
evaluation,  or  12.6  sec  per  function/derivative  evaluation  (23  function  evalu¬ 
ations).  The  integration  has  been  performed  with  250  integration  steps,  and 
it  may  be  possible  to  reduce  the  time  by  using  forward  differences.  However, 
the  convergence  characteristics  would  deteriorate. 

6. 2  Plane  Change  Problem 

For  this  problem,  a  good  initial  guess  of  the  final  time,  the  angle  of 
attack,  the  bank  angle,  and  the  ignition  time  have  been  found  to  be  t^.  -  2733  see  , 
n  -  40  deg  ,  u  *  60  deg  ,  and  t^  •  1624  deg  .  The  corresponding  trajectory 
is  shown  in  Table  10,  from  which  it  is  seen  that  the  orbital  inclination  is 
if  ■  22.2  deg  and  the  constraint  residuals  are 
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C  =  hf/ 608,000  -  1  - 

C2  =  Vf/23,500  -  1  -  .882155E-2 

C3  =  Yf  -  .109798E-1  (63) 

c4  =  y7  -  -o 

C5  :  y8  “  *° 

The  last  two  constraints  are  the  inequality  constraint  on  the  angle  of  attack 
(a  <_  40  deg)  and  an  inequality  constraint  on  the  bank  angle  (p  0) .  In  the 
solution  of  the  minimization  problem,  the  bank  angle  at  the  final  point  tried 
to  go  to  a  large  negative  value.  The  last  inequality  constraint  has  been 
included  to  keep  this  from  happening. 

The  set-up  for  the  plane  change  problem  (IPROB  -  2)  is  similar  to  that 
of  the  reentry  problem.  An  additional  parameter,  the  ignition  time,  is 
present  making  a  total  of  twelve  (N  -  12).  Here,  there  are  five  constraints 
(M  -5),  of  which  three  are  equality  constraints  (K  ■  3).  At  first,  the  values  of 
EPS  had  been  chosen  to  be  l.E-4  *  X(I)  .  However,  because  of  the  extreme 
sensitivity  of  the  problem  to  the  ignition  time,  it  has  beer,  necessary  to  use 
EPS (I)  *  l.E-5  *  X(I)  .  Finally,  the  constraint  scale  factors  C(M  +  I)  for 
the  altitude,  the  velocity,  the  flight  path  inclination,  and  the  ineouality 
constraints  have  been  set  equal  to  .009  ,  .009  ,  .011  ,  .001  ,  and  001  , 
respectively.  The  first  three  values  are  rounded  constraint  residuals, and  the 
last  two  have  been  chosen  to  produce  relatively  small  weights. 

A  summary  of  the  iteration  process  is  presented  in  Tables  11  through  14,  and 
the  converged  trajectory  is  illustrated  in  Figures  4  through  6.  After  four 
cycles,  the  constraints  are  satisfied  to  four  significant  figures.  The 
converged  angle  of  attack  history  stays  near  40  deg  (maximum  C  )  at  the 
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Table  12.  Maximum  Plane  Change:  Cycle  2  Results 
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Figure  5.  Maximum  Plane  Change:  0  ,  $  ,  and  h  versus 


(deg ) 


Figure  6.  Maximum  Plane  Change:  V  ,  ,  and  1  versus 


high  altitudes  ond  <j<  creases  to  13  deg  (maximum  C  /C_)  at  the  lov  altitudes. 

L  U 

On  the  other  hand,  the  bank  angle  is  near  70  deg  during  the  entry  part  of 
the  trajectory  and  decreases  to  0  deg  during  the  exit  part.  The  maximum  plane 
change  is  i^  »  33.7  deg  ,  and  the  load  factor  does  not  exceed  0.5  g's. 

The  CPU  time  for  the  four  cycles  is  4175  sec  on  the  CYBER  170/750 
computer.  This  amounts  to  approximately  0.5  sec  per  function  evaluation  and 
10  sec  per  derivative  evaluation. 

For  this  problem,  the  Lagrange  multipliers  v  do  not  seem  to  be  approach¬ 
ing  a  limit.  A  possible  reason  is  that  the  problem  is  so  sensitive  to  the 

value  of  t.  that  additional  accuracy  is  needed  in  the  numerical  Integration 
b 

to  achieve  better  convergence  accuracy.  Also,  in  generating  numerical  deri¬ 
vatives,  the  same  perturbation  size  is  used  for  each  variable.  More  consistent 
results  can  be  achieved  by  tailoring  the  perturbation  size  to  each  variable. 

The  behavior  of  the  bank  angle  of  the  final  point  can  be  explained  by 
assuming  that  the  earth  is  not  rotating.  The  Initial  point  of  the  trajectory 
is  located  on  the  equator,  and  the  vehicle  is  trying  to  put  the  velocity  vector 
on  the  orbital  plane  with  the  highest  inclination.  As  long  as  the  longitude, 

6  , is  less  than  180  deg,  the  bank  angle  for  maximizing  the  inclination  is 
positive.  However,  once  the  vehicle  crosses  the  equator  into  the  southern 
hemisphere,  it  must  bank  in  the  opposite  direction  (p<0)  to  obtain  additional 
inclination.  By  restricting  the  bank  angle  to  positive  values,  the  maximum 
inclination  should  be  obtained  when  0^  ■  180  deg  and  p  -  0  .  The  optimal 
trajectory  presented  here  has  these  features,  as  can  be  seen  from  Table  14. 

On  the  other  hand,  it  should  be  possible  to  increase  the  Inclination  by 
relaxing  the  inequality  constraint  to  u  ^  -90  deg  . 
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SECTION  VII 


DISCUSSION  AND  CONCLUSIONS 

The  maximum  crossrange  and  the  maximum  plane  change  problems  have  been 
formulated  as  parametar  optimization  or  nonlinear  programming  problems.  An 
existing  code  for  solving  the  nonlinear  programming  problem  vith  the  augmented- 
Lagranglan  method  has  been  selected  and  modified  slightly  for  solving  these 
problems.  The  derivatives  needed  for  the  optimization  method  have  been  computed 
numerically.  As  a  consequence,  the  entire  program  can  be  viewed  as  a  model 
for  solving  virtually  any  trajectory  optimization  problem. 

A  difficulty  in  solving  trajectory  optimization  problems  lies  in  the 
model  used  for  the  problem.  For  example,  in  performing  the  one-dimensional 
search,  it  is  possible  to  change  the  parameters  so  much  that  resulting  trajec¬ 
tories  go  to  extremely  high  altitudes  or  extremely  low  altitudes.  This  has 
caused  an  exponential  overflow  in  the  atmosphere  subroutine.  In  another  case, 
the  resulting  trsjectorles  ended  up  in  a  vertical  dive,  causing  the  numerical 
integration  to  blow  up  because  cf  the  singularity  in  the  heading  angle  equation. 
Initially,  this  difficulty  had  beer  eliminated  by  removing  the  singularity, 
which  means  increasing  the  number  of  differential  equations  by  one  and  Increasing 
the  numerical  complexity  of  the  program  (time  becomes  a  dependent  variable). 
Later,  by  additional  modification  of  the  optimization  code,  it  has  become 
possible  to  monitor  the  flight  path  inclination  and  terminate  the  one-dimensional 
search  when  y  gets  too  small.  Another  difficulty  arose  by  trying  to  modify 
the  drag  polar  so  that  the  inequality  constraint  on  the  angle  of  attack  could 
be  eliminated.  The  optimization  process  has  found  a  way  to  use  the  change 
to  its  advantage  and  has  produced  an  unrealistic  trajectory.  Finally,  in  the 


46 


plane  change  problem,  the  bank  angle  at  the  final  pi  nt  wants  to  go  to  a 
large  negative  value.  At  this  point,  it  is  not  understood  what  aspect 
in  the  model  makes  it  useful  for  this  to  happen.  Hence,  an  inequality 
constraint  has  been  imposed  on  the  bank  angle  to  keep  it  nonnegative. 

Once  the  modeling  problems  have  been  eliminated,  the  reentry  problem 
is  relatively  easy  to  solve.  This  is  not  the  case  for  the  plane  change 
problem,  as  can  be  seen  from  the  number  of  iterations  and  the  computer 
time  required.  Possible  reasons  include  the  number  of  constraints  Involved, 
sensitivity  of  the  problem  to  ignition  time,  and  numerical  accuracy,  either 
in  the  Integration  or  in  the  computation  of  derivatives.  With  regard  to 
constraints,  the  plane  change  problem  is  easy  to  solve  if  the  bank  angle  is  the 
only  control  and  the  altitude  is  the  only  constraint  at  the  final  point.  Adding 
the  other  final  conditions  reduces  greatly  the  progress  achieved  on  every 
iteration.  Also,  including  the  angle  of  attack  as  a  control  reoulres  the  use 
of  the  inequality  constraint  on  a  . 

The  particular  method  for  computing  the  derivatives,  central  differ¬ 
ences,  has  been  chosen  because  it  is  easy  to  modify  the  model.  (If  the 
derivatives  are  computed  from  differential  equations,  any  change  in  the 
model  produces  a  change  in  the  derivative  equations.)  Although  it  has  not 
been  done,  it  is  possible  to  cut  the  computer  time  roughly  in  half  by  only 
integrating  over  that  part  of  the  trajectory  which  is  affected  by  a  per¬ 
turbation.  Also,  the  computer  time  can  be  reduced  further  by  using  forward 
differences,  or  one  function  evaluation  per  derivative  calculation.  This 
reduction  in  accuracy  will  affect  the  accuracy  with  which  the  constraint 
can  be  satisfied. 
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The  nanner  In  which  the  optimization  problems  have  been  solved  la 
sufficient  to  demonstrate  the  ability  of  the  optimization  code  to  handle 
Much  problems.  Questions  concerning  the  number  of  nodes  needed  to  obtain 
an  accurate  solution  and  the  placement  of  these  nodes  are  still  open,  but 
they  can  be  answered  with  a  moderate  amount  of  computer  time. 


48 


Mbs 


REFERENCES 


1.  Bryson,  A.  E. ,  Jr.,  and  Ho,  Y-C,  Applied  Optimal  Control,  John  Wiley 

and  Sons,  New  York,  1975. 

2.  Fletcher,  R.,  "An  Ideal  Penalty  Function  for  Constrained  Optimization", 

J.  Inst.  Maths.  Appl..  Vol.  15,  pp.  319-342,  1975. 


49 


APPENDIX  A 

DOCUMENTATION  ASSOCIATED  WITH 
THE  OPTIMIZATION  CODE 
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Marvell  Subroutine  Library 


VF01  A/AD 


1.  Purpose 

To  find  the  minimum  of  s  function  F(*)  of  n  variables  x  subject  to  both 
equality  constraints  c^(_x)  •  0  i  ■  i,2,..7.k  and  inequality- constraints 
Ci<x)  >0  1  •  k  t  1 , •  •  • , m  Ik  •  0  or  k  •  m  is  allowed  but  k  must  be  4  n). 
Derivatives  of  all  func‘.  Ions  with  respect  to  x  must  be  provided,  both  the 

vector  (dF/0x^ ,  dF/dx2 . ,3F/<)xn)  and  the  matrix  whose  ith  column  is 

(dej/dx-j,  dcj/dx2,..'.tdcj/dxn)  for  1  a  i,2..,,m.  These  functions  and 
derivatives  must  be  specified  in  a  user  subroutine  VF01B  (see  section  4). 

An  initial  estimate  of  the  solution  must  be  provided  which  need  not  be 
feasible.  The  subroutine  allows  advantage  to  be  taken  of  the  possibility  that 
some  of  the  constraints  are  linear,  and  also  of  certain  other  types  of  infor¬ 
mation  about  the  problem,  if  available.  If  all  the  constraints  are  linear, 
the  use  of  VFOiA  is  not  most  efficient,  and  one  of  the  LA  or  VE  routines 
should  be  used.  The  method  is  a  penalty  function  -  Lagrangian  method  (see 
section  8)  and  VFOIA  calls  VA09A  to  carry  out  the  associated  unconstrained 
minimisations. 

2.  Argument  List 

CALL  VF01AIN, M.K.X, EPS, AKMIN.DFN.KAXFN.IPRl, IPR2, IW, MODE) 

N  An  INTEGER  set  to  the  number  of  variables  n  (N  *2). 

M  An  INTEGER  set  to  the  total  number  of  constraints  m  (M 

K  An  INTEGER  set  to  the  total  number  of  equality  constraints  k. 

X  A  ,i4_aL  array  of  N  elements  in  which  the  initial  estimate  of  the  solution 

must  be  set.  VFOiA  returns  the  solution  x  in  X. 

EPS  A  REAL  array  of  N  elements,  in  which  the  tolerances  for  the  unconstrained 
minimizations  must  be  set.  EPS (I)  should  be  set  so  that  EPS(I)/X(I)  v 
AKMIN,  roughly  speaking. 

AKKIN  A  REAL  number  in  which  the  relative  error  tolerance  required  in  the 
constraint  residuals  must  be  set.  VFOIA  will  exit  when  max|lc^(x)  I  / 
scaling  factor  for  cjJ  {  AKMIN  for  the  active  constraints  jij  (see 
section  7). 

PFN  A  REAL  number  in  which  the  likely  reduction  in  F(x)  must  be  set.  This 
is  done  in  the  same  way  as  for  VA09A,  -  see  the  VA09A  specification 

sheet. 

MAXFN  An  INTEGER  in  which  the  maximum  number  of  calls  of  VFQ1B  on  any  one 
unconstrained  minimization  must  be  set.  Roughly  speaking  2  or  3  times 
MAXFN  calls  of  VF01B  are  likely  to  be  made  altogether. 

IPR1  An  INTEGER  controlling  the  frequency  of  printing  from  VFOiA.  Printing 

occurs  every  XPRl  iterations,  except  for  details  of  increases  to  the  0^ 
which  are  always  printed.  No  printing  at  all  occurs  (except  for  error 
diagnostics )  if  IPR1  •  0.  All  printing  controlled  by  IPR1  ia  suitably 

annotated.  VFOIA  1 
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IPR2  An  INTEGER  controll  lng  the  frequency  of  printing  from  VAO'.V..  U>ti2 
should  be  aet  as  described  in  the  VA09A  specification  sheet. 

IW  An  INTEGER  giving  the  amount  of  storage  available  In  COMMON/ v/rOlL/W( . ). 
Sat  to  2'j00  unless  wishing  to  change  the  restrictions  (see  Section  5). 

MCOE  An  INTEGER  controlling  the  mode  of  operation  of  VF01A.  If  any  positive 
definite  estimate  is  available  of  the  heaaian  matrix  of  the  penalty 
function,  set  |M00e|  -  2  or  3,  otherwise  eet  |moue|  «i  (see  VA09A 
specification  sheet).  If  estimates  of  the  and  0^  parameters  are 
available  (aee  section  8)  aet  MODE  <  0,  otherwise  set  MOOE  >0.  A 
normal  setting  for  e  one-off  j6b  with  no  information  available  is 
MtOE  -  l. 

3.  The  named  COMMON  areas 


Certain  named  COMMON  area*  must  be  declared  and  aet  on  entry  to  VFOlA. 


COMtON/VFOlC/CdSO) 


COf*ON/Vf'OiF/GC(25, 50 ) 


COMMON/ VTO  1G/T(  150 ) 


COMMON/ vro l I/G2P< 325 ) 


Set  scale  factors  <>0)  for  the  constraints  in 
C(1),C(2),....,C(M).  Choose  the  magnitude  of 
these  scale  factors  to  qlve  an  Indication  of  the 
magnitude  of  the  constraints  evaluated  about  the  Initial 
approximation  x.  If  any  constraints  are 
violated  by  an  amount  graatar  in  modulus  than 
that  which  is  set,  then  the  setting  is  increased 
accordingly,  n vast  seals  factors  are  transferred 
to  C(M.l),  C(M.2),.....,C(2M)  by  VFOlA. 

Set  the  derivatives  of  any  linear  constraints  on 
entry  rather  than  in  VF01B.  This  is  the  moat 
efficient  and  the  numbers  are  not  disturbed.  The 
manner  of  setting  is  described  in  section  4. 

If  MCOE  <  0  is  used,  then  set  the  oarameters 
e1(e2,  ....  dK  in  T<1),  T(2),...,T(M)  and  the 
parameters  ^.cr-,  ...,<rn  in  T(M*1),T(M»2) , ... ,T(2M). 

The  meaning  of  these  parameters  may  be  found  in 
the  report  TP552  -  aee  section  8. 

If  I MOOE I  «  2  or  3  set  the  estimated  heaaian 
matrix  of  the  penalty  function  in  G2P(1>,..., 
G2P(N*(N«l)/2).  The  manner  of  setting  is  that 
described  in  the  specification  sheet  of  VA09A  under 
the  heading  MOOE. 


4.  The  user  subroutine  VF01B 

The  user  must  declare  a  subroutine  headed 

SUBROUTINE  VF01B(N,M,X) 

REAL  X(l) 

CCMMON/VF01C/F 
COMMON/ VP01D/G( 50) 

COMMON/ VF0l£/C( 150 ) 

COMMON/VF01F/GC<  25 , 50 ) 

Thia  subroutine  must  taXe  x,  the  vector  supplied  in  X(1),...,X(N)  and  set 

F(x)  in  f'l  C.(x),...  >^(£7  in  C(l) . C(M>1  WF/dx.  ,  ...<3F/dxn)x  in  G(l),... 

G(N)j  and  aet  7dc^/dx^,...eJci/JXn)x  in  GC(1,I) , . . .  ,GC( N, I )  for  aTl  I  ■  1,2, •• 
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(Excepting  the  linear  constraints  wMun  st.mld  i*.  r,et  on  entry,  as  the 
numbers  are  constant].  Some  time  m  ■■■•  ,<lso  ho  s  .v—i  If  required  by  also  Inclu¬ 
ding  COMMON/ Vf01G/T(  150)  and  by  not  -’valua'  ii.g  -  ,n  ,i ) . GC(M,I )  when 

c ( I >  J,  T ( I )  for  any  I  >  K.  Note  thm.  the  opr-lmm-  values  F(x),  (d  F/d*lt . . . , 
dt/dxn,  etc.  are  left  in  these  named  Common  aro-i:.  on  exit  from  VF01A.  Note 
also  that  In  the  double  precision  version  tuo  usoi  subroutine  name  is  VP01BD 
and  there  is  a  0  appended  to  the  named  COMMON  areas. 

5.  Redefining  named  COMMON  areas 

Local  storage  for  VF01A  is  through  named  COMMON  areas.  These  have  been 
set  on  the  assumption  that  k  4  ?i  and  M  ^  50.  If  it  is  desired  to  r«nove 
either  or  both  of  these  restrictions,  then  it  is  necessary  to  increase  the 
storage  available  in  some  or  all  of  these  areas.  This  can  be  done  by  defining 
the  named  COMMON  areaa  in  the  users  MAIN  with  the  increased  storage  settings, 
in  which  case  the  extra  storage  will  be  effective  throughout  the  whole  program. 
The  complete  list  of  named  COMMON  used  by  VF01A  and  the  corresponding  values 
of  N  and  M  are  at  follows. 


COMMON/ VF01C/F,M,K, IS, MK,NU 
D/G ( 50 ) 

E/C (150) 

F/GC(2S,  50)* 
G/T(1S0) 

H/GP(50) 

I(G2P(325) 

J/V(50) 

K/WW(1S0) 
L/W(2S00 ) •• 
M/ZZ(100) 
N/LT(100) 


Independent  of  N  and  M 
2N 
3M 


N,M 

3M 

p  (p  •  max(M.N)) 

N*(N*l)/2 

P 


?P 

?M 


Notes:  ^ 

•  On  Increasing  M,  when  N  <  25,  redefine  GC  with  bounds  (N.M)  not  (25, M). 

VF01A  has  been  coded  under  this  assumption,  as  it  requires  less  storage. 

(VF01A  treats  GC  as  a  single  suffix  array). 

••For  M  very  large,  storage  locations  may  be  prohibitively  large.  However 
it  Is  very  unlikely  that  this  amount  of  storage  will  actually  be  needed  by 
VF01A  (no  problem  has  yet  been  encountered  for  which  more  than  (?N)2  locations 
have  been  required).  Hence  in  these  circumstances,  either  declare  w  with 
(2N)2  locations  (or  whatever  can  be  spared),  and  set  the  lnteqer  IW  to  this 
number  in  the  calling  sequence  of  VF01A.  If  more  locations  are  required, 
then  VF01A  will  stop  and  print  out  the  storage  required. 

6.  General 

yat  aL-cmaw 

'toSUtViSii 


Other  routines: 


Inout/Cutput : 


Restrictions: 

VF01A  3 


named  COMMON  only  -  see  section  5. 

5K  words  unless  N  or  M  is  redefined,  when  it  is  not 
usually  more  than~4$N2  ♦  nm  *  0(max  (M,N))  words. 

Calls  VF01B  (user  subroutine),  Vf\)lZ  (private),  VA09A, 
MC11A,  MCUE.  VA09A  calls  MCUB  in  addition. 

No  input,  all  output  on  stream  6  (line  printer),  con¬ 
trolled  by  user  througn  IPR1  and  2PR2. 

N  $  25  and  M  $  50  but  can  be  lifted  -  see  section  5. 
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Pate  of  routine! 


August  1973 


7.  Accuracy 

It  may  be  that  VF01A  Is  unable  to  achieve  the  accuracy  requested  In  the 
parameter  AKMIN.  In  this  case  a  diagnostic  Is  printed.  To  find  the  cause 
of  this,  first  examine  tna  print  out  of  the  VA09A  iteration.  If  this  Is 
anomalous  (pd  A->  0  for  Instance)  suspect  a  mistake  in  the  programming  of  VF01B 
particularly  in  obtaining  derivatives.  If  VA09A  seems  O.K. ,  then  other 
possible  causes  are  (1)  there  Is  no  feasible  point  (in  which  case  oA  ->  >  and 
ci  *  const  A  0),  (11)  EPS  has  been  set  too  large  relative  to  AKMIN,  (ill) 

AKMIN  has  been  set  to  sstall  relative  to  the  machine  precision,  (iv)  the 
problem  is  too  ill-conditioned. 

8.  Method 

That  described  in  R.  Fletcher.  «An  ideal  penalty  function  for  constrained 
optimization",  C.S.S.2,  December,  1973.  The  penalty  function  for 
inequalities  1$ 


d(x,j6,ff)  «  F(jt)  »  ^  £  o-j  (c^lx)  “  Qi)2 

and  each  iteration  involves  minimizing  6(x)  for  fixed  6_,o.  After  each  itera¬ 
tion  the  8  and  £  parameters  are  varied  so  that  the  sequence  of  minima  £a|<r 
tends  to  The  solution  of  constrained  problem.  The  value  of  the  producT^ 

©lOj  tends  to  the  1th  Lagrange  multiplier  of  the  problem.  Any  information 
about  Lagrange  multipliers,  or  about  the  hessian  of  6  can  usefully  be 
incorporated. 

Convergence  is  guaranteed  (in  exact  arithmetic)  and  this  implementation 
of  the  method  can  be  expected  to  converge  at  a  second  order  rate. 


October,  1973 
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Harwell  Subroutine  library 


V  AO  9  A/A  D 


t .  Purpose 

lb  find  the  minimum  of  a  function  F(x)  of  several  variables,  given 
that  the  gradient  vector  (3F/9x),0F/3x2,.T.,0F/^xn)  can  be  calculated. 

The  subroutine  replaces  VAOIA  to  which  it  is  superior  in  various  ways  (see 
section  5),  and  should  be  used  whenever  derivatives  can  be  evaluated  readily. 
It  should  however  not  be  used  either  if  storage  space  is  at  a  premium  (use 
VA08A)  or  if  the  function  is  a  sum  of  squares  (use  VA07A).  The  subroutine' 
complements  VAOfiA,  the  latter  requires  four  times  the  storage,  and  some 
comparisons  (R.  Fletcher,  A.E.R.E.  Report,  R712S  (1972))  indicate  that 
VAOat  is  marginally  slower  and  more  affected  by  round  off  error.  As  VAOOt 
is  more  difficult  to  use,  it  is  suggested  that  VA09A  should  be 
used  in  the  first  instance  on  ar\y  problem.  If  VA09A  falls  then  VAO^  should 
be  tried  as  it  is  guaranteed  to  converge  if  the  effect  of  rounding  errors  can 
be  neglected. 

2.  Argument  List 

GALL  VA09A(FUNCT,N,X,F,G, H,W, DFK, EPS, MODE, MAXFN,IPRINT,IEXIT) 

FUNCT  An  IDENTIFIER  of  the  users  subroutine  -  see  section  3. 

N  An  INTEGER  to  be  set  to  the  number  or  variables  (N  >2). 

X  A  REAL  ARRAY  of  N  elements  in  which  the  current  estimate  of  the 

solution  is  stored.  An  initial  approximation  must  be  set  in  X 
on  entry  to  VA09A  and  the  best  estimate  obtained  will  be 
returned  on  exit. 

F  A  REAL  number  in  which  the  best  vulue  of  F(x)  corresponding  to 

X  above  will  be  returnee. 

G  A  REAL  ARRAY  of  N  eianents  in  which  the  gradient  vector 

corresponding  to  X  above  will  be  returned.  Not  to  be  set  on 
entry. 

H  A  REAL  ARRAY  of  N*(N+i)/2  elements  in  which  an  estimate  of  the 

2 

hessian  matrix  3  F/(3x.3x  )  is  stored.  The  matrix  is  represented 

in  the  product  form  LDL  where  L  is  a  lower  triangular  matrix  with 
unit  diagonals  and  D  is  a  diagonal  matrix.  The  lower  triangle 
of  L  Is  stored  by  columns  in  H  excepting  that  the  unit  diagonal 
elements  are  replaced  by  the  corresponding  elements  of  D.  The 
setting  of  H  on  entry  is  controlled  by  the  parameter  RODE  (q.v.). 

W  A  REAL  ARMY  of  3*N  elements  used  as  working  space. 
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DFN  A  REAL  number  which  mat!  In'  -et  so  as  co  give  VA09A  an 

estimate  of  the  likely  reuui  Uon  to  be  obtained  in  F(x). 

DFN  is  used  only  on  the  fii~t  iteration  so  an  order  o? 
magnitude  estimate  will  s,ullice,  The  information  can  be 
provided  in  different  ways  depending  upon  Uic  sign  of  DFN  which 
should  be  set  in  one  of  the  following  vmys: 

DFN>0  the  setting  of  DfN  itself  will  he  taken  as 
the  likely  reduction  to  be  obtained  in  F(x). 

DFN=0  it  will  be  assumed  that  an  estimate  of  the 

minimum  value  of  F(j$)  has  been  set  in  argument 
F,  and  the  likely  reduction  in  F(x)  will  be 
computed  according  to  the  ini tial”fu ration  value. 

DFN<0  a  oultiple  | DFN |  of  the  modulus  of  the  initial 
function  value  will  be  taken  as  an  estimate  of 
the  likely  reduction. 

EPS  A  REAL  ARRAY  of  N  elements  to  be  set  on  entry  to  the  accuracy 

required  in  each  elonent  of  X. 

VODE  An  IOTEGBt  which  controls  the  setting  of  the  initial  estinte 
of  the  hessian  matrix  in  the  parameter  H.  The  following 
settings  of  MDDE  are  permitted. 

MDCE=1  An  estimate  corresponding  to  a  unit  matrix 
is  set  in  H  by  VA09A. 

MODE-; 2  VA09V  assumes  that  the  hesBian  matrix  itseir 
has  been  set  in  H  by  columns  of  its  lower 

T 

triangle,  and  the  conversion  to  LDL  form 
Js  carried  out  by  VAO&A.  The  hessian  matrix 
must  be  positive  definite. 

MODELS  VA09A  assumes  that  the  hessian  matrix  has  been 

set  in  H  in  prodict  form.  This  is  convenient  when 
using  the  H  matrix  from  one  problem  as  an  initial 
estimate  for  another,  in  Wiich  case  the  contents 
of  H  are  passed  on  unchanged. 

M\XFN  An  INTEGER  set  to  the  maximum  number  of  calls  of  FUNCT  permitted. 

IPRINT  An  INTEGER  controlling  printing.  Printing  occurs  every 
| IPRINT)  iterations  and  also  on  exit,  in  the  form 

Iteration  No,  No  of  calls  of  FUNCT, IEXIT  (on  exit  only) 

Function  value 

X( 1) ,X(2) , . . . ,X(N)  8  to  a  line  (5  in  VA08AD) 

G(1),G(2),...,G(N)  8  to  a  line  (5  in  VA08AD) 

The  values  of  X  and  G  can  be  suppressed  on  interne diate 
iterations  by  setting  IPRINTcO.  All  intermediate  printing  can 
be  suppressed  by  setting  IPRINT=M\XFN+"i .  Ail  printing  can  be 
suppressed  by  setting  IPRINT=0, 
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I  EXIT  An  INTEGER  giving  the  r  -ion  for  exit  from  VA09A.  Thi* 
will  be  let  by  VA09A  as  .  ollowg 

IEXIT=0  (mode=2  only).  The  estimate  of  the  hessian 
matrix  is  not  positive  definite. 

IEXIT=I  The  normal  exit  in  v>  ich  |QX(1)|<EPS(I)  for 
aii  1=1, 2,... N,  where  DX(I)  is  the  change  in 
X  on  an  iteration. 

T 

IEXIT=>2  C  □ X^D.  Not  possible  without  rounding  error. 
Probable  cause  is  that  EPS  is  set  too  miall 
for  computer  word  length. 

IEXIT=3  FUNCT  called  MAXFN  times. 

3.  User  Subroutine 


The  user  must  provide  a  subroutine  headed 
SUBROUTINE  XXX(N,X,F,C) 

REAL  X( 1 ) ,G( 1 )  (REAL’S  in  VAOSAD) 


where  XXX  is  an  ioentifier  chosen  by  the  user. 

This  subroutine  should  use  the  variables  x  supplied  in  X(l), 
X(2),.,.,X(N)  to  evaluate  the  function  and  gradient  vector  and  place 
them  in  F  and  G(1),G(2) , . . . ,G(N)  respectively.  .  XXX  must  be  passed  to 
VA09A  as  VA09A' s  first  argument,  see  section  2,  and  appear  in  an 
EXTERNAL  statesient  in  the  program  that  calls  VA09A. 

4.  General 

Use  of  CO*CN  :  none 

Workspace:  N*(N>l)/2  words  +  4N  words  provided  by  the  user  in 

H  and  W. 


Other  routines:  calls  MC11A,  6,  E. 

Input/Output:  controlled  by  the  user  through  IPRINT.  Ail 

output  is  on  stream  6  (line  printer). 

Restrictions:  none 


System  dependence:  none 

Date  of  routine:  April,  1972, 
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5.  Method 


The  method  used  is  a  qua  si. -flew  ton  method  described  by  Fletcher 
(Computer  Journal,  Vol.  13,  p,317,  1970),  and  is  a  modification  of 
earlier  methods  of  this  tyr  such  as  that  implemented  by  VAOIA.  The 
method  is  superior  to  that  ■  »  VAOIA  on  three  counts, 

(1)  It  uses  a  formula  to  update  the  hessian  approximation  H 
which  has  proved  to  be  more  efficient  and  reliable. 

(2)  It  uses  a  'crude'  line  search  which  has  been  shown  to  be 
more  efficient  than  an  'accurate'  line  search. 

T 

(3)  It  represents  H  by  the  product  LDL  ,  which  enables  the 
positive  definiteness  of  H  to  be  guaranteed,  even  in  the 
presence  of  round-off  error. 


Harwell  Subroutine  Library 


MCI  I  A/AD 


i 

! 


t ,  Purpose 

MCI1A  ia  a  subroutine  for  use  in  problems  which  involve  the  addition 
or  subtraction  of  rank  one  matrices  c  zzT  to  positive  definite-or  seal- 
definite  symmetric  matrices  A  stored  In  factored  fora  A  =  LDL,  such  that 
the  resulting  N  x  N  matrix 

A  =‘A  *  <J  tt 


is  also  known  to  be  positive  definite  or  seal-definite.  Note  that  L  la 
lower  triangular  with  4^*1,  and  D  la  diagonal  with  d^  >,  0.  Apart  froa 

Its  obvious  use  in  updating  Matrices  which  reaein  strictly  positive 
definite,  MC11A  can  be  used  (1)  to  acauaiiate  a  sun  of  rank  one  teras  into 
an  initial  matrix  A=0,  (11)  to  carry  out  projection  end  allied 

operation  on  A  which  redice  the  rank,  and  (ill)  to  update  Matrices  A  of 
rank  k  <  n  where  it  is  known  from  other  considerations  that  the  rank 
reaains  unchanged.  All  these  operations  preserve  the  correct  rank  and  are 
not  seriously  affected  by  round-off  error,  Ihe  Method  is  that  described 
by  M.J.D.  Powell  and  R.  Fletcher  (1973),  ’On  the  updating  of  LDLT 
factorizations' ,  T.P.  519. 

The  Matrix  A  is  represented  using  the  niniaal  storage  of  N*(N+l)/2 
elements  where  N  is  the  distension  of  the  problem.  To  facilitate 
operating  with  A,  a  number  of  independent  subroutines  have  been  provided, 
written  as  entry  points  MC11B/C/. ../F.  These  operations  include  reducing 
a  matrix  to  its  factors,  sultiplylng  out  the  factors,  operating  with  the 
factors  of  A  on  a  vector  t  to  obtain  either  Az  or  A~'z,  and  replacing  the 
factors  of  A  by  the  matrix  A-*.  These  facilities  are  described  in  more 
detail  in  section  4. 

2* 

CALL  UMIACA.N.Z.SIO.W.IR.MC.EPS) 

A  A  REAL  one  dimensional  array  of  N*(N+l)/2  elements  in  vtiich  the  matrix 
T 

A=LDl  mst  be  given  in  factored  fora.  The  order  in  which  el  an  en  tv  of 

L  and  D  are  stored  is  <Ve2t*e3l**>*’Wd2,e»'**’'<kl2' . 

n_i The  factors  of  the  matrix 

•v  T 

A  =  A  az  j  will  overwrite  those  of  A  on  exit. 

N  An  XNIGtXR  (N>l)  which  must  be  set  to  the  dimension  of  the  problem, 

Z  A  REAL  one  dimensional  array  of  N  elements  in  which  the  vector  z  must 

be  set.  The  errey  Z  Is  overwritten  by  the  routine. 

8IQ  A  REAL  variable  in  which  the  scalar  c  sust  be  set.  SIO  is  not 
restricted  to  11.,  but  if  SICkO  then  it  aust  be  known  from  other 
considerations  that  X  it  positive  definite  or  semi -definite,  apart 
from  the  effects  of  round-off  error. 
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W  A  REAL  Array  of  N  element*.  If  SIG>0  then  W  is  not  used,  and  the 
name  of  any  one  dimensional  array  can  be  inserted  in  the  calling 
sequence.  If  SIG<0  then  W  is  used  as  work  space.  In  addition  for 
SIG<0  it  nay  be  possible  to  save  time  by  setting  in  W  the  vector 
£  defined  by  Lv=£.  The  ways  in  which  this  can  occur  are  described 
under  MC  belowT 


IR  An  INTCCQl  to  be  set  so  that  |IR|  is  the  rank  of  A.  If  the  rank  of  A 
is  ejected  to  be  different  from  that  o£  A,  set  IR4O.  On  exit  from 
MCI  1A,  IR(j.O)  will  contain  the  rank  of  A. 


IK  An  INTEGER  to  be  set  only  when  SIG<0,  as  follows.  If  the  vector 
v  defined  by  has  not  been  calculated  previously,  set  MM), 

Tf  MCI  IE  lus  been  used  previously  to  cslculate  A-1z,  then  visa  by¬ 
product  or  thic  calculation  and  is  stored  in  the  ^parameter  of  MCI  IE. 
In  this  case  transfer  v  to  the  W  parameter  of  MCI  1A  and  aet  MCst. 

If  l  tus  been  calculated  as  z  *  Au  for  some  arbitrary  vector  u  using 
MCiTd,  then  again  v  Is  a  by-prodict  of  the  calculation  and  is” 
available  in  the  W~peraoeter  of  MCI  ID.  In  this  case  (or  any  other 
in  which  v  ia  known)  set  v  in  the  W  parameter  of  MCltA  and  set  MC*2. 

EPS  A  REAL  variable  to  be  aet  only  when  SIC<0  and  A  is  expected  to  have 
the  same  rank  as  A.  ^  In  ggrt^in  ill-conditioned  cases  a  non-zero 
diagonal  element  of  D  (X=LDLr)  might  become  so  small  as  to  be 
Indeterminate,  Two  courses  of  action  are  possible.  One  is  to 
introduce  a  small  perturbation  in  order  that  X  keeps  tlie  same  rank 
as  A.  This  ia  the  normal  course  of  action  and  is  achieved  by  setting 
EPS  equal  to  the  relative  machine  precision  e.  e  is  -  I0~7  in  single 
length  arithmetic  and  -  lO"1®  in  double  length  on  the  IBM  370.  The 
other  course  or  action  is  to  let  the  rank  of  X  be  one  less  than  the 
rank  of  A.  This  is  achieved  by  setting  EPS  equal  zero. 

3.  Oenerai 

Use  of  OCSaCN :  none 


Workspace: 


N*<N-*-l  )/2  ♦  N  *  N  words  provided  by  the  user  in  A,Z  and  W. 
If  SIG>0  W  need  not  be  supplied. 


Entry  points: 


MCI  IB/C/.  ../IF  -  see  section  4, 


Other  routines:  none. 


Inout/OutiHJt: 

Timing: 


none. 

2 

One  call  of  MC11A  requires  *■  n  multiplications,  unless 

2 

SIGcO  and  M(aO  when  the  figure  is  *“  iJj  n  .  One  call  of 

3  / 

any  of  MCI  IB, C  or  F  requires  ~  n  /6  suitiplicaticns.  One 

2 

call  of  either  MCI  ID  or  E  requires  -  n  /2  sultiplications. 


Restrictions:  none. 

System  dependence:  none. 

Pete  of  routine:  Jamary, 


'1 

■3 

1973.  j 

) 


4.  Outer  Entry  point* 


1 


Other  entry  points  are  provided  to  facilitate  operating  with  A 
which  is  stored  in  compact  form.  In  all  of  these  A  is  a  REAL  one 
dimensional  array  of  N*(N-fi),^  elements,  and  N^l  is  an  integer  set  to  the 
dimension  of  the  problan,  Each  entry  point  is  essentially  an  independent 
sv ^routine,  and  could  be  taken  out  and  written  as  such  if  desired* 

MCHB  -  factorize  a  positive  definite  syuaetric  matrix  given  in  A. 

CALL  MCI  1B(A,N,IR) 

A  Mist  contain  the  elesients  of  A  in  the  order  a, a,.,, a  „a  ,,... 

•m . VllWV,'W  ° 

successive  columns  of  its  lower  triangle).  On  exit  A  will  be  over¬ 
written  by  the  factors  L  and  D  in  the  form  described  in  section  2  under  A. 

N  Order  of  the  matrix  A. 

IR  An  INTEGER  set  by  MC1IB  to  the  rank  of  the  factorization.  If 

the  factorization  has  been  performed  successfully  IR=N  will  be  set. 

If  IR<N  then  the  factorization  has  failed  because  A  is  not  positive 
definite  (possibly  die  to  round-off  error).  In  this  case  the 
factors  of  a  positive  a  end-definite  matrix  of  rank  IR  will  be  found  In 
A.  However  the  results  of  this  calculation  are  unpredictable,  and 
MCHB  should  not  be  used  in  an  attempt  to  factorize  positive  rani- 
definite  matrices. 

T 

MC11C  -  miltiply  out  the  rectors  LDL  to  obtain  A. 

CALL  MCIIC(A.N) 

A  Mist  be  set  in  the  factored  form  described  in  section  2  under  A. 

On  exit  the  factors  will  be  overwritten  by  the  explicit  matrix  A, 
the  order  of  the  elements  being  that  described  for  input  to  MCI  IB. 

N  Order  of  the  matrix  A. 

MCI  ID  -  calculate  the  vector  z*=Az  where  A  is  In  factored  form. 

CALL  MCI  !D(A,N,Z,W) 

A  Mist  be  set  In  factored  form. 

N  Order  or  the  matrix  A. 

Z  A  REAL  array  of  N  elements  to  be  set  to  the  vector  z.  On  exit 
z  contains  the  vector  i** Az. 

W  A  REAL  array  of  N  elements  which  is  set  by  MCI  ID  to  the  vector  v 

defined  by  Lv=z*.  If  this  vector  is  not  of  interest,  replace  W  by 
Z  In  the  calling  sequence  to  obviate  the  need  to  supply  extra  storage. 
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fcCI  It:  -  calculate  the  vector  where  A  la  in  factored  fora. 

CALL  MCI 1E(A,N(Z,M.  IR) 

A  Mist  be  set  in  factored  form. 

N  Order  of  the  natrix  A. 

Z  A  REAL  array  of  N  elements  to  be  set  to  the  vector  j.  On  exit  Z 
contains  the  vector  z*=A_1  z, 

W  A  REAL  array  of  N  elements  widen  is  set  by  IC11E  to  be  vector  v 

defined  by  Lv»j.  If  this  vector  is  not  of  interest!  replace  W  by  2 
in  the  calling  sequence  to  obviate  the  need  to  supply  extra  storage. 

1R  An  INTEGER  which  cust  be  set  to  the  rank  of  A. 

1C11F  -  calculates  the  e>q>lici  t  matrix  A-1  fro*  the  factors  of  A. 

CALL  IC11F(A,N,IA) 

A  Hist  be  set  in  factored  form.  ICI1F  will  overwrite  this  by  the 

elements  of  the  inverse  matrix  A**1,  in  the  order  Br}>A2t,**',aNN  “ 
for  MCI  IB. 

N  Order  of  the  matrix  A. 

IR  An  INTEGER  which  oust  be  set  to  the  rank  of  A. 

Notes.  (i)  SCI  IF  should  not  be  used  to  solve  equations,  in  which  case 
MCI  IE  should  be  used.  MCllF  is  intended  for  applications  in  which  the 
explicit  elanents  of  A"’  must  be  examined,  for  example  in  the  use  of 
variance-covariance  matrices.  (ii)  MCI  IE  and  F  will  RETUSN  without 
calculation  unless  UUN. 


A 

1 


January,  1975 
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VE04A/AD 


I ,  Purpose 


To  find  the  value  of  % 


xi 

*fi) 


which  minimizes  a  '-.aadratic  function 


QOg)  of  n  variables  subject  to  upper  and  lower  bounds  on  some  or  all  of 
the  variables.  Q(g)  is  defined  by 

n  n  n 

Q(«)  •  W**  -  "  *5  1  Y  XiAiJXj  ■  ^\*i 

ivl  J=l  1=1 


where  A  ia  aynmetric  (A^sA^),  and  the  bounjs  are  of  the  form  I^ax^au^.  It 
is  peruiissible  to  let  i^=u,  if  required,  ancl  A  is  not  restricted  to  being 
positive  definite.  Tne  subroutine  calculates  the  solution  .»  =  the  mininum 
value  Q(£) ,  and  the  gradient  g(£)  (note  £(x)=Ax-£).  This  problem  is  u 
special  case  of  quadratic  programming  for  which  a  subroutine  VE02A  exists. 
VE04A  is  more  efficient  and  reliable  for  sol'  Ing  problems  of  this  special  rorm. 

An  application  of  the  subroutine  is  to  (weighted)  linear  least  squares 
data  fitting  subject  to  bounds.  If  it  is  required  to  minimize 


S(2t)  =  (B*-y)TW(B*-y) 


(2) 


subject  to  the  above  bounds,  where  E  is  an  mxn  matrix  m>n  and  W  is  an  mxm 

diagonal  rrotrix  of  weights  Wu>0,  inpn  set  A=2BTWB  ard  .£  =  2BTP[y.  Statistical 

calculations  for  this  problem  are  described  in  section  5,  including  an  additional 

entry  point  to  enable  the  variance-covariance  iratrix  to  be  calculated. 

2.  The  Anainent  list 
♦ 

CALL  VE04A  (N,A,  IA,B  ,BL,0U  ,X,Q,  LI' , K,G) 

N  an  ItfTEGER  which  must  be  set  by  the  user  to  the  nunber  of  variables. 

A  a  REAL,  two  dlnereioml  urrqy  ,  each  dincnsion  at  least  N;  the 

elements  in  the  upper  triangle  A(l,J)  T*JtN  must  be  set  by  the 

V30l*A  1 


63 


H.lifSwm-***4' 


user  to  the  cor  reloading  A^  in  (I},  end  will  remain  untouched 
ty  the  aubroutine.  Elerents  A(I,J)  N»I>J  arc  used  as  working  apace. 

1A  an  INTEGER  giving  the  first  dimension  of  A  in  the  statement  which 
assigns  space  to  A. 

B  a  REAL  array  of  at  least  N  elements.  The  user  must  set  B(I)  to  the 
bA  in  (l).  B  is  not  overwritten  by  VEOtA. 

BL  a  REAL  array  of  at  least  N  elements.  The  user  oust  set  BL(I)  to 

the  lower  bound  ly  on  the  Ith  variable.  If  the  bound  is  non-existent, 
set  l y  to  s  very  small  nuaber  like  -1E75.  BL  is  not  overwritten  by 
VEOtA. 

BU  a  REAL  array  of  at  least  N  elenents.  The  user  must  set  Bu(I)  to  the 
upper  bound  u^  on  the  1^*  variable.  If  the  bound  ia  non-existent, 
set  Uj  to  a  very  large  number.  BU  Is  not  overwritten  by  VE04A. 

X  a  REAL  array  of  at  least  N  elements.  VE04A  returns  the  solution 
iy  inX(l). 

Q  a  REAL  variable  in  which  VEOtA  returns  the  value  of  Q(£). 

LT  an  INTEGER  array  of  at  least  N  elements,  set  by  VE04A  to  a 
permutation  of  the  integers  1,2,. ..,N  (see  K  and  C  below) 

K  an  INTEGER  set  by  VE04A  to  u»e  nunfcer  of  free  variables  at  the 

solution  £  (those  not  on  Umir  bounds).  These  are  the  variables 
LT( I ) ,  LT(2),...,LT(K). 

G  A  REAL  array  of  at  least  3*N  elements,  G(1 ),.... ,G(N)  are  set  by 
VE04A  to  the  gradient  g(£). .  C  is  indirectly  addressed  so  that 
C(I)  contains  the  gradient  with  respect  to  the  LT(I)  variable, 
whence  0(1 COO  will  be  found  to  be  zero.  G(N*1),..., 

G(3*N)  are  used  by  VE04A  as  working  *>ace. 


2  VdOU 


6  4 


3. 


Geirral  Information 


Use  of  CCMrtON:  none 

2 

Workspace:  approx  words  (fair  of  A)  +  2n  words  in  G  and  n  integers 
in  LT. 

i her  routines:  none 
Entry  points:  VE04B-see  section  5 
Input/Output:  none 
System  dependence:  none 

Timing.  The  time  required  depends  upon  how  many  free  variables  k  there 
are  at  the  solution.  Typical  f.  ures  for  (k/n,  nunt>er  of 
multiplications)  are  (.1,  n3/l2),  (.5,n3/6),  (,75,n3/2). 
Restrictions:  none 
Date  of  Routine:  April  1973. 

4.  Me thod 


That  of  Fletcher  R.  and  Jackson,  M.l*.,  (1973),  "Minimization  of  a  quadratic 
function  of  many  variables  subject  only  to  lower  and  upper  bounds",  T. P.528. 

This  method  conbines  generality  (any  A),  efficiency  (times  comparable  to  those 

7 

required  for  a  factorization  of  A)  and  stability  (uses  partial  LDL  factorizations). 
5.  Statistical  Calculatii~-ns 

When  a  sum  of  squares  is  being  minimized  as  in  (2),  then  certain  statistical 

quantities  can  readily  be  calculate-!  Firstly,  of  course,  S(§)  is  given  by 
T 

Q(^)  .  _y  W(y,  If  it  is  assumed  that  Ur  bound  variables  located  by  VE04A  are 
exact  in  the  underlying  model,  then  an. estimate  of  the  residual  variance  Is 
given  by  S(^)/(m-k) .  To  estimate  variatxcs  and  covariances,  an  additional 
entry  point  is  provided.  This  calculates  { A }  1  where  1 A |  indicates  the  submatrix 
jAAj  1  where  i  and  j  index  only  the  free  variables.  Tlx:  appropriate  variance- 
covariance  matrix  for  the  free  variables  is  then  o*  1 A i ”  1 .  Estimates  of 
standard  deviations  of  the  free  variables  are  given  by  the  square  roots  of  the 
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diagonal  clinunts  of  this  matrix.  Because  the  bound  variables  are  known 
exactly,  Uuy  have  sero  variance  and  covariance. 

The  entry  point  VE048  is  essentially  written  as  a  separate  subroutine. 
It  calculates  1a}*"1  and  is  used  as  follows: 

CALL  VE04B  (N.A.IA.C.K) 

N,A,IA,G  and  K  must  be  passed  on  unchanged  after  exit  from  VE04A. 
VE048  sets  the  following: 

A  The  off-diagonal  elements  of  |a]-1  are  set  in  A(I,J)  for 

J<lclC.  The  elements  are  indirectly  addressed  so  that  A(I,J) 

contains  {a}”*  where  r=LT(I)  and  s»LT{J). 

r  ,a 

C  The  diagonal  elements  of  IaI”'  are  set  in  G(N+1) . .  C(N+K). 

They  are  again  indirectly  addressed  so  thntG(N*I)  contains 
-1 

Ulr>r  where  r=LT(I). 


2nd  May  1973 
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SUBROUTINE  VFO 1 A( N , M , K , X , EPS , AKMIN , DFN . MAXFN , IPR 1 , IPR2 , IW , MODE) 
REAL  X( 1 ) ,EPS( 1 ) 

COMMON/ VFO 1C/F, MM, KL. IS, HK.NU 
COMM0N/VF01D/GC50) 

COMMON/ VF0 1 E/ C( 1 50  > 

COMMON/VFO 1 F/GC( 1 250 ) 

C0MM0N/VF01G/TO50) 

COMMON/VFO 1H/GP(50) 

COMMON/VFO 1I/G2P( 325) 

COMMON/ VFO 1J/V( 50) 

COMMON/VFO 1K/WW( 150) 

COMMON/VFO 1L/W( 2500) 

COMMON/VFO 1M/ZZ( 100) 

COMMON/VFO 1N/LT( 100) 

COMMON/TIME/ TMAX , TO , T 1 , IRS , MINS , AK , ITN , IFN 

COMMON/INOM/INOM 

EXTERNAL  VF01Z 

DATA  BOUND/ 1.E+8/ 

1000  FORMATOOI1)) 

1001  F0RMAT(8E15.7) 

NU=MAXO(25.N) 

IF(M.GT.50)NU=N 

IX=N 

ICS=M 

ICB=M+M 

IS=M 

IL=IS+M 

IP=M 

ILT=M 

NN=N*(N+1 )/2 

MM=M 

KL=K 

Rsl. 

MK=0 

INOMsI 

CALL  VFOIB(N.M.X) 

DFsDFN 

IF( DFN . LT.OEO ) DF=ABS( DFN*F) 

IF( DFN.EO.OEO)DFsF 
IF(DF.LE.0.)DF=1. 

IF ( MODE . LT . 0 ) G0T05 

DO  2  Is  1  ,M 
CC=C(I) 

IF(I.GT.K)CC=AMIN1(CC,0E0) 

IF(  ABS(  CC)  .GT.C(  ICS«-I)  )C(  ICS+I)sABS(CC) 

2  CONTINUE 
DO  3  1*1, M 

T( IS*I)=2E0#DF/C(  ICS+I)*#2 

3  T( I ) sO . 

5  CONTINUE 
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IF(  IPR1 .  EQ.0)GOTO4 
IF(MODCMINS.IPR1).NE.0)G0TO4 
PRINT  1002 

1002  FORMATC' OENTRY  TO  VF01A'///'0C0NSTRAINT  SCALE  PARAMETERS  ARE') 
PRINT  1001 ,(C<ICS*I), 1=1 ,M) 

4  CONTINUE 
MD=IABS( MODE) 

8  CONTINUE 

00  9  1=1, NN 

9  W( I)=G2P( I) 

IF( IPR1 . EQ.0)G0T07 
IF(M0D(MINS,IPR1).NE.0)G0T07 
PRINT  1003, MINS 

1003  F0RMAT( //// ' OOUTER  ITERATION  NUMBER  IS1, 13) 

PRINT  1004 

1004  FORMATC 'OX(I) ' ) 

PRINT  1001, (X(I), 1=1, N) 

PRINT  1005 

1005  FORMATC 'OTHETA( I)') 

PRINT  1001 ,<T(I) ,1=1 ,M) 

PRINT  1006 

1006  FORMATC OSIGMAC I)') 

PRINT  1001  .(TUS+I),  1=1,  M) 

7  CONTINUE 

IFCIRS.EQ.1)  GO  TO  82 
ITN=0 
IFN=0 
82  IRS=0 

CALL  VA09AC VF012.N.X, PHI ,GP,W,WW, DF, EPS.MD.MAXFN, IPR2, IEXIT) 

IFCTI-TO.LT. THAX)  GO  TO  80 

DFN=DF 

DO  81  1=1, NN 
8i  G2PCI)=W(I) 

RETURN 
80  CONTINUE 
INOM= 1 

CALL  VF01BCN.M.X) 

MD=3 

AKK=0. 

DO  10  1=1, M 

IFC I . GT . K . AND . CC I ) .EQ.O.)  C(I)=1.E-8 
CC=CCI) 

IFC I.GT.K. AND.CC I) .GE.TC I) ) CC=AMIN1 ( CC.OEO) 

TCI)  =  TCI)»TC  IS+I) 

WW(  I)=ABS(  CC)/CC  ICS-*- 1) 

IFCWWCI) . GT . AKK ) AKK= WWC I ) 

10  CONTINUE 

IFC IPR1 . EQ.O)GOTOl6 

IFC MODC MINS, IPR1) ,NE.O)GOTOl6 

PRINT  1007 
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1007  FORMATCOEXIT  FROM  VA09A ' / ' OLAGRANGE  MULTIPLIERS  USED  IK  VA09A' ) 
PRINT  1001, (T(I), 1=1, M) 

PRINT  1008 

1008  FORMAT* 'OLARGEST  SCALED  CONSTRAINT  VIOLATION’/ 

1’  THIS  ITERATION,  BEST  ITERATION') 

PRINT  1001,AKK,AK 
PRINT  1009 

1009  FORMATCOCONSTRAINT  RESIDUALS') 

PRINT  1001, (C(I), 1=1, M) 

PRINT  1010 

1010  FORMATC ' OSCALED  CONSTRAINT  VIOLATIONS') 

PHINT  1001 ,(WW(I) ,1=1 ,M) 

16  CONTINUE 

IF( IEX IT . EQ . 0 . OR . IEXIT . EQ . 3 ) GOT020 
IF(  AKK .  LE . AKMIN)GOT020 
IF( AKK.GE.AK)G0T01 1 
DO  15  1=1 ,  NN 
15  G2P(I)=W(I) 

DO  17  1=1, M 

IF ( I.GT.K. AND. T( I) . EQ.OEO . AND. C( I) .GT.0E0)G0T017 
ZZ(IP+I)=-T(IS+I)»C(I) 

IF(  I.GT.K. AND. ZZC IP+I) .LT.-TC I) )ZZ(IP*I)=-T(I) 

17  CONTINUE 
IF(MINS.EQ.I) G0TO40 
GOTO 18 

11  CONTINUE 

DO  14  1=1, M 

IF{  WW(  I ) . LE . AK . OR . C( ICB* I ) . GE . 4E0»WW( I ) ) G0T01 4 
DS=9E0*T( IS+I) 

T(IS+I)=1E1*T(  IS+I) 

IF( IPR1 .NE.O) PRINT  101 1 ,I,T<IS*I) 

1011  FORMATC 'OSIGMAC .13,')  INCREASED  TO  '  ,E15.7) 

DO  12  J=1,N 

12  V(J)=GC((I-1 )*NU+J) 

CALL  MC11 A(G2P,N,V,DS,V,N,N,DS) 

14  CONTINUE 

18  CONTINUE 

DO  13  1=1, N 

IF  C ABSCXC I)-G( IX+I) ) .GT.EPS( I) )G0T021 

13  CONTINUE 
PRINT  1013 

1013  FORMATCO REQUESTED  ACCURACY  NOT  OBTAINED') 

20  CONTINUE 

DO  2012  1=1 , NN 
2012  G2P( I)=W( I) 

IF( IEXIT. EQ.0)PRINT  2000 

2000  FORMATC OMATRIX  SET  IN  G2P  BY  USER  IS  NOT  POSITIVE  DEFINITE') 

IF( IPR1 . EQ.O) RETURN 
PRINT  1012 

1012  FORMATCOBEST  SOLUTION  OBTAINED'/'  -,-,(G( I) ,  1=  1  ,N) ' ) 
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PRINT  1001, F.(G(I), 1=1. H) 

RETURN 
21  CONTINUE 

IF(  AKK .  LT .  AK )  G0T040 
DO  32  1*1, H 
32  V(I)*T(IUI> 

GOTO70 
NO  CONTINUE 
MKtO 
KK=0 

DO  41  1=1 ,M 
T(IUI)=T(I) 

C(ICB4.I)=WW(I) 

IF(  I.GT.K.AND.T(IL+I) . EQ.OEO. AND.C(  I) .GT.OEO)GOT041 

KKeKK+1 

LT( ILT+KK)=I 

GP(KK)=-BOUND 

IF(I.GT.K)GP<KK)=-T<IUI) 

V(KK)= BOUND 
ZZ(KK)=-C(  I) 

41  CONTINUE 

IF(KK . EQ.0)GOT020 
DO  42  1=1  ,N 

42  G( IX+I)=X( I) 

KKK=KK* (KK*1 )/2 
II=MAXO(KKK+NN,KK#KK) 

IF(II.LE.IW)G0TO50 
PRINT  2001,11 

2001  FORMAT! 'OINCREASE  STORAGE  IN  COMMON/VFOIL  TO' .17 ELEMENTS' ) 
RETURN 

50  CONTINUE 
II=IW-KKK 

DO  53  1=1 . KK 
LI=LT( ILT+I) 

DO  51  JJ=1,N 

51  X(JJ)=GC( (LI-1 )*NU*JJ) 

CALL  MCI 1E(W,N,X,DUM1 ,X,N, IDUH.DUH2) 

DO  53  J«1.I 
LJ*LT( ILT+J) 

Lz  0. 

DO  52  JJ=1,N 

52  Z«Z«-X(JJ)»GC(<LJ-1)#NU*JJ) 

ZIsIZ^I 

53  W(II).Z 
JJsIW-KKK 
II«0 

DO  56  1*1, KK 
DO  55  J=1,I 
JJ.JJ^I 

55  N(  II*J)*W(  J«!) 
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56  II=II+KK 

CALL  VE04A(KK,W,KK,ZZ,GP,V,T,Z,LT,JJ.W) 

IF( IPR1 .EQ.C)G0T059 
IF<M0D(MlNS,iPR1) .NE.0)G0T059 
PRINT  1020, KK 

1020  F0RMATU4,'  ACTIVE  CONSTRAINTS,  NUMBERED*  > 

PRINT  1000,(LT( ILT+I) ,1=1 ,KK) 

PRINT  1021 

1021  FORMATS OLAGRANGE  MULTIPLIER  CORRECTIONS  FOR  ACTIVE  CONSTRAINTS') 
PRINT  1001, (TCI), 1=1, KK) 

59  CONTINUE 

DO  60  1  =  1  ,M 

60  V( I)  =  T( IL+I) 

DO  62  1=1, KK 
LI=LT(  ILT+I) 

V(LI)=V(LI)+T(I) 

Z=4E0*ABS( (T( I)-ZZ( IP+LI) )/ZZ( IP+LI)) 

IF(Z.LE.1E0)G0TO62 
DS=( Z-1£0)*T( IS+LI) 

T( IS+LI )=Z*T( IS+LI) 

IF( IPR1 . NE.O) PRINT  101 1 , LI, T( IS+LI) 

DO  61  J=1  ,N 

61  GP(J)=GC((LI-1)»NU+J) 

CALL  MCI 1A(G2P,N,GP,DS,GP,N,N,DS) 

62  CONTINUE 
AK=AKK 

70  CONTINUE 

DO  71  1=1 ,M 

71  T(I)=V(I)/T(IS+I) 

DO  72  1=1,  N 

72  X( I)=G( IX+I) 

DF=1E50 

MINSsMINS+1 

G0T08 

END 
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SUBROUTINE  VA09A(FUNCT, N.X.F.G.H, W.DFN.EPS, MODE, MAXFN.IPRINT, 
1  IEXIT) 

REAL  X( 1 ) ,GC 1 ) ,H( 1 ) ,W( 1 ) ,EPS( 1 ) 

COMMON/T IME/ TMAX , TO , T 1 , I RS . MINS , AK , ITN , IF* 

COMMON/ IGAMMA/ IGAMMA 
IF( IPRINT.NE.O)PRINT  1000 
1000  FORMATC//' ENTRY  TO  VA09AV) 

NNsN«(Nf1)/2 

IG=N 

IGGsN+N 

ISsIGG 

IEXITsO 

IRsN 

DFrDFN 

I F  ( MODE . EQ . 3 ) GOT0 1 5 
IF(  MODE. EQ.2)G0T010 
IJ=NN«-1 
DO  5  1=1, N 
DO  6  Jsl.I 
IJ=IJ-1 
6  H(IJ)=0. 

5  H(I J)  =  1 . 

GOTO 15 
10  CONTINUE 

CALL  MC1 1  B(H,N,  D1 ,  D2,D3,  IK,  ID1  ,D4) 

IF(IR.LT.N) RETURN 
15  CONTINUE 
Z=F 

CALL  FUNCT(N.X.F.G) 

IFN=IFN*1 

IF(DFN.EQ.O.)DF=F-Z 
IF( DFN . LT .0 . ) DF=ABS< DF»F) 

IF( DF.LE.0.)DF=1 . 

20  CONTINUE 
DFN=DF 

CALL  SECOND(Tl) 

IF(T1-T0.GE.TMAX)G0  TO  90 
IFCIPRINT.EQ.0)G0T021 
IF  ( M0D( ITN , IPRINT) . K£ ,0)G0T021 
PRINT  1001, ITN, IFN 

1001  F0RMAT(2**I5) 

PRINT  1002, F 

1002  FORMAT((8E15.7)) 

IF{ IPRINT.LT. 0/G0TO21 
PRINT  1002, (X(I) ,1*1 ,N) 

PRINT  1002, (G(I>, 1*1, ») 

21  CONTINUE 
lTN*XTHr1 
DO  22  I«1 ,N 

22  WUG*I’;»G<I) 
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CALL  MC1 1 6(H,N,G,D1 ,V, IR, ID1 ,  D2) 

os=o. 

DO  29  r=1  ,N 
W( IS+I)=-G(I> 

29  GS=GS-G{ I)*W( IG+I) 

IEXITs2 

IF(GS.GE.0.)G0T092 

GSOeGS 

ALPHAs-2 . *DF/GS 

IF( ALPHA. GT. 1.0)  ALPHAsl.O 

DF=F 

TOTsO. 

30  CONTINUE 
IEXITs3 

IF< IFN . EQ . MAXFN)G0TO92 

ICONsO 

IEXITd 

DO  31  1=1  .  N 

Z=ALPHA*W( IS-*- 1) 

IF(ABS(Z) .GE.EPS(I)) ICONs 1 

31  X(I)=X(I)«.Z 

CALL  FUNCTlN.X.FY.G) 

IFNsIFN+1 

IFC IGAMMA. EQ.O)  GO  TO  33 
DO  34  Isl.N 
ZsALPHA*W(IS«-I> 

34  X(I)sX(I)-Z 
ALPHA=0.1* ALPHA 
PRINT  35. ALPHA 

35  FORMAT(1X.»ALPHAs*,F10.5) 

GO  TO  30 

33  CONTINUE 
GYSsO. 

DO  32  >,H 

32  GYS=GYS*G(I)*WvIS+I) 
IF(FY.GE.F)G0T040 

IF( ABS(GYS/GSO) .LE. ,9)G0T050 
IF( GYS . GT . 0 . ) G0TO40 
TOTsTOT+ALPHA 
Z*10. 

IF(GS. LT.GYS) Z.GYS/<  GS-UYS) 

IF( Z.GT . 10 . ) Zs 10 . 

ALPHAs ALPHA* Z 
FcFY 
GSsGYS 
GOT030 

40  CONTINUE 

DO  41  Isl.N 

41  X(  I)=X<  I)-ALPHA*W(  IS+t) 
IF(IC0N.EQ.0)G0T092 
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Z: 3 • * ( F-F Y ) / ALPHA+GYS+GS 
2ZsSQRT( Z**2-GS*GYS) 

Zsl .-(GYS+ZZ-Z)/(2.*ZZ+GYS-CS) 

ALPHA: ALPHA*Z 

GOT030 

50  CONTINUE 
ALPHAbTOT*.  ALPHA 
FaFY 

IF(ICON.EQ.0)GOTO90 
DF=DF-F 
DGS=GYS-GSO 
DO  51  Isl.N 
W(  IGG>I)=G(I) 

51  G(I)=-W( IG+I) 

IF(  DGS+ALPHA*GSO . GT . 0 . ) GOTO60 
C  COMPLEMENTARY  DFP  FORMULA 
SIG=1 ,/GSO 
IRs-IR 

CALL  MCI  1A(H,N,G,SIG,Vf,IR,  1,0.) 

DO  52  Isl.N 

52  G( I)sW( IGG+I)-W(  IG4-I) 

SIGsl ./( ALPHA*DGS) 

IRs-IR 

CALL  MCI 1A( H.N.G.SIG.W, IR.0,0.) 
G0T070 

60  CONTINUE 
C  DFP  FORMULA 

ZZ=ALPHA/( DGS-ALPHA*GSO> 

SIGs-ZZ 

CALL  MCI 1A(H,N.G.SIG.W,IR,1, IE-13) 
Z=DGS»ZZ-1 . 

DO  61  Isl.N 

61  G(I)sW( IGG+I)*Z*W(  IG+I) 

SIG=1 ./(ZZ«DGSB»2) 

CALL  MCI  1A(H,N,GtSIG,Vf, IR.0,0.) 

70  CONTINUE 

DO  71  Isl.N 

71  G(I)sW(IGG4.I) 

GOTO20 

92  CONTINUE 
DO  91  Isl.N 
91  G(I)«W(IGd) 

90  CONTINUE 

IF(IPRINT.EQ.O)RETURN 
PRINT  1001 , ITN, IFN, IEXIT 
PRINT  1002, F 
PRINT  1002, (X(I), Isl.N) 

PRINT  1002, (0(1), I.i.N) 

RETURN 

END 
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SUBROUTINE  VF01Z(H,X,PHI,GPHI) 
REAL  X( 1 ) ,GPHI( 1 ) 

COMMON/VF01C/F,M,K, IS.MK.NU 
COMMON/VF01D/G<50) 

COMMON/ VF01E/C( 150) 
COMM0N/VF01F/GC( 1250) 
COMMON/VF01G/T(150) 

COMMON/ IGAMMA/ IGAMMA 
IF(MK.EQ.1)CALL  VF01B  (N,M,X) 
MK=1 

IF( IGAMMA. EQ.1)  RETURN 
PHI=0. 

DO  10  1=1, N 

10  GPHI( I)sG(I) 

DO  12  1=1  ,M 
CC=C( I)-T(I) 

IF( I ,GT.K)CC=AMIN1 ( CC.OEO) 

Y=T( IS+I)*CC 
IF( Y. EQ.0E0)GOT012 
PHI=PHI+Y»CC 
DO  11  J=1,N 

1 1  GPHI( J)=GPHI( J)+Y*GC( ( 1-1 )*NU+J) 

12  CONTINUE 
PHI=.5E0»PHI+F 
RETURN 

END 
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SUBROUTINE  MC 11 A( A , N . Z , SIG . V, IR , MK , EPS) 
DIMENSION  A( 1 ) ,Z( 1 ) ,W< 1 ) 

C  UPDATE  FACTORS  GIVEN  IN  A  BY  SIG»Z«ZTRANSPOSE 
IF(N.GT.1)G0T01 
IR=  1 

A(1)=A( 1)fSIG  *Z( 1 )**2 
IF( A(1). GT.O.) RETURN 
A( 1 )=0. 

IR=0 
RETURN 
1  CONTINUE 
NPsN+1 

IF( SIG. GT.O. )GOTO40 

IF(  SIG. EQ.O. .  OR .  IR .  EQ .  0 ) RETURN 

TI=1 ./SIG 

IJsl 

IF( MK. EQ.0)GOTO10 
DO  7  1=1  ,N 

IFCA(IJ) . NE.O.)TI=TI-fW( I)*,2/A(IJ) 

7  IJ=IJ-*-NP-I 
G0TO20 

10  CONTINUE 

DO  11  1=1 , N 

11  W(I)=Z(I) 

DO  15  1=1. N 

IP=I+1 

V=W(I) 

IF( A( IJ) .GT.O. ) GOTO 12 
W(I)=0. 

I J=I J+NP— I 
GOTO 15 

12  CONTINUE 
TI=TI*V*<>2/A(IJ) 

IF( I . EQ. N)G0T014 
DO  13  J=IP,N 
IJ=IJ+1 

13  W(J)=W(J)-V»A(IJ) 

14  IJ=IJ+1 

15  CONTINUE 

20  CONTINUE 
IFUR.LE.O  )GOT021 
IF(TI .GT.O. )G0TO22 
IF(MK-1>40,40,23 

21  TIsO. 

IR=-IR-1 

G0T023 

22  TI=EPS/SIG 

IF( EPS. EQ.O.) IRaIR-1 

23  CONTINUE 
MM=1 


77 


MC1 1 A  PAGE  2 


TIMsTI 
DO  30  1=1,  N 
J=NP-I 
IJ=IJ-I 

IF(  A( I J ) . NE .0 . ) TIK*TI-W<  J)**2/A( IJ) 

W(J)=TI 
30  TIsTIM 
G0T041 
AO  CONTINUE 
MM=0 

TIM=1  ./SIC- 
A1  CONTINUE 
IJsl 

DO  66  1=1 , N 

IPsI+1 

V=Z(I) 

IF( A( IJ) .GT.0.)G0T053 

IF( IR.GT.O  .OR.SIG.LT.O. .OR.V.EQ.O. )G0T052 
IR=1-IR 

A( IJ)=V**2/TIM 
IFC I. EQ.N) RETURN 
DO  51  J=IP,N 
IJ=IJ*1 

51  A( IJ)=Z( J)/V 
RETURN 

52  CONTINUE 
TI=TIM 
IJ=IJ+NP— I 
G0T066 

53  CONTINUE 
AL=V/A(IJ) 

IF(MM)5A.5A,55 

5A  TI=TIM+V»AL 
G0T056 

55  TI=W( I) 

56  CONTINUE 
R=TI/TIM 
A(IJ)=A(IJ)*R 
IF(R.EQ.0.)G0T070 
IF(I.EQ.N)GOT070 
B*AL/TI 

IF(R.GT.A.)G0T062 
DO  61  J=IP,N 
IJ=IJ+1 

Z<J)=Z(J)-V«A(IJ) 

61  A(IJ)  =  A(IJ)+B»Z(J) 

GGT06A 

62  GH.TIH/TI 

DO  63  J«IP,N 
IJ»IJ«-1 
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Y=A(IJ) 

A(  IJ)sB*Z(«J)+Y*GM 

63  Z(J;=Z(J)-V*y 

64  CONTINUE 
TIMsTI 
IJsIJ*1 

66  CONTINUE 
70  CONTINUE 

IF(IR.LT.O)IR=-IR 

RETURN 

C  FACTORIZE  A  MATRIX  GIVEN  IN  A 
ENTRY  MC1 1 B 
IR=N 

IF(N.GT.1)G0TO100 
IF(A(1).GT.0.)RETURN 
A( 1 )sO  . 

IRsO 

RETURN 

100  CONTINUE 
NPsN+1 
11*1 

DO  10^4  I=2,N 
AArA(II) 

NIs II+NP-I 

IF( AA.GT.O .) GOTO 101 

A( II)sO . 

IRrIR-1 

IIsNI+1 

GOTO104 

101  CONTINUE 
IPsII+1 
II=NI* 1 
JK=II 

DO  103  IJ=IP,NI 
V=A(IJ)/AA 
DO  102  IK=IJ,NI 
A( JK)aA(JK)-A(IK)»V 

102  JK=JK-*.1 

103  A(IJ)»V 

104  CONTINUE 
IF(A(II).GT.0.)RETURN 
A(II)aO. 

IR*IR-1 

RETURN 

C  MULTIPLY  OUT  THE  FACTORS  GIVEN  IN  A 
ENTRY  MCI  1C 
IF(N.EO.I) RETURN 
NPsN+1 
II*N*Np/2 
DO  202  NIPs2,N 
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JKsII 

MIaII-1 

II=II-NIP 

AAsA(II) 

IP=II«-1 

IF<  AA . OT . 0 . ) G0T0203 
DO  204  IJsIP.NI 
204  A(IJ)aO. 

GOTO202 
203  CONTINUE 

DO  201  IJ=IP,NI 
V=A( IJ)*AA 
DO  200  IK=IJ,NI 
A(JK)  =  A(JK)«-A(IK)»V 

200  JK= JK+1 

201  A( I J)*V 

202  CONTINUE 
RETURN 

C  MULTIPLY  A  VECTOR  Z  BY  THE  FACTORS  GIVEN  IN  A 
ENTRY  MCI  ID 
IF(N.GT.1)GOT0300 
Z(1)=Z(1)*A<1) 

W(1)=Z(1) 

RETURN 

300  CONTINUE 
NP=N+1 
11=1 
N1=N-1 

DO  303  1=1. N1 
Y=Z( I) 

IF( A( II) . EO.O . )GOTO302 

IJsII 

IP*I+1 

DO  301  J=IP,N 
IJaIJ+1 

301  Y=Y+Z(J)»A(IJ) 

302  Z( I)=Y#A( II) 

W(I)=Z(I) 

303  II=II+NP-I 

Z(  N)  =  Z(  N)  •'AC  II) 

W(N)*Z(N) 

DO  311  Kal.NI 
IaN-K 

Ilell-NP+I 

IF( Z( I) . EQ.O. )G0T031 1 
IP»L»1 
I  Jail 
YsZ(I) 

DO  310  J*IP,N 
IJaIJ+1 
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310  Z(  J )=Z( J)+A( IJ)*2( I) 

311  CONTINUE 
RETURN 

C  MULTIPLY  A  VECTOR  Z  BY  THE  INVERSE  OF  THE  FACTORS  GIVEN  IN  A 
ENTRY  MCI  IE 
IF(IR.LT.N) RETURN 
WC 1 )=2< 1 ) 

IF(N.GT.1)GOT0400 
Z(  1 )  sZ(  1  )/A(  1 ) 

RETURN 

400  CONTINUE 

DO  402  1=2, N 
IJ:1 
11=1-1 
V=Z(I) 

DO  401  Jsl.II 
V=V-A( I J)*Z( J) 

1101  IJsIJ+N-J 
W(I )=V 
402  Z(I)=V 

Z(N)=ZCN)/A(IJ) 

NPsN-f  1 

DO  Nil  NIP=2,N 
I=NP-NIP 
II=IJ-NIP 
V=Z( I)/A( II) 

IP=I*1 

IJ=II 

DO  410  JsIP.N 
IIsIUI 

410  V=V-A(II)»Z(J) 

411  Z(I)=V 
RETURN 

C  COMPUTE  THE  INVERSE  MATRIX  FROM  FACTORS  GIVEN  IN  A 
ENTRY  MCI  IF 
IF( IR. LT.N) RETURN 
A(  1  )  =  1 ./ A( 1 ) 

IF(N.EQ.I) RETURN 

NP=N+1 

NlsN-1 

11=2 

DO  511  I»2.N 
A( II)*-A( II) 

IJ=II«.1 

IF(I.EQ.N)G0TO502 

DO  501  J=I,N1 

IK*  II 

JK«IJ 

VsA(IJ) 

DO  500  K»I,J 
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JK=JK+NP-K 

V=V*A<IK)*A(JK) 

500  IK=IK+1 
A( IJ)s~V 

501  IJsIJ+l 

502  CONTINUE 
A(IJ)=1 ./A(IJ) 

AA=A(IJ) 

IJsI 

IP=I«-1 

NI=N-I 

DO  511  J=2,I 

V=A(IJ)»AA 

IK=IJ 

KsIJ-IPW 

I1=IJ-1 

NIPsNI+IJ 

DO  510  JK=K,I1 

A(JK)=A(JK)+V*A(IK) 

510  IK=IK-fNIP-JK 
A( I J)=V 

511  IJ=IJ«.NP-J 
RETURN 
END 
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SUBROUTINE  VEONA(  N , A, IA, B, BL, BU.X ,Q,LT,K,G) 

DIMENSION  A(  IA,  1 )  ,B(  1 ) ,BL( 1 ) ,&U( 1 ) ,X( 1 ) ,LT( 1 ) ,G(  1 ) 

IS=N 

IASsN 

IVsN 

ICAC=N+N 
IDsICAC 
DO  9  1=1.  N 

9  G(I)s-B(I) 

DO  10  1=1, N 
X(I)=0. 

LT  ( I )  =  I 

G(  ICAC+I )=A(  1,1) 

IF(0. .GE.BL( I) .AND.O. .LE. BU( I) )GOTO10 
I F  ( 0 . .  LT .  BL(  I ) )  X  ( I )  =  BL(  I ) 

IF(0.  .GT.BIK  I)  )X(  I)=BU(  I) 

DO  12  J=1 ,1 

12  G( J  )=G( J)+A( J,I)*X(I) 

IF( I.EQ.NJGOTOIO 
11=1+1 

DO  11  JsII.N 
11  G(J)=G(J)+A(I.J)#X(I) 

10  CONTINUE 
K=0 

K1  =  1 

20  CONTINUE 
IOUTsO 
DEL=0. 

DO  21  I=K1,N 
LI=LT(I) 

IF(  X( LI) . EQ. BL(LI) .AND.G(I) ,GE.0.)GOTO21 
IF(X(LI) .EQ.BU(LI) .AND.G(I) .LE.O. )G0T021 
IF( G( I) . LT.O. )G0T022 
Z=X(LI)-BL(LI) 

J=1 

G0T023 

22  CONTINUE 
Z=B'J(LI)-X(LI) 

J=0 

23  CONTINUE 

IF( G( ICAC+I) .LE.0.)GOTO2N 
BETAsABS(G(I) )/G(  ICAO I ) 

IF(BETA.GE.Z)G0T024 

ZsBETA 

Ds.5#Z»ABS(G(I)) 

Js-1 
G0TO26 
2«  CONTINUE 

DsZ*( ABS(G( I) )-.5*Z*G( ICAC+I) ) 

26  CONTINUE 
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IF( 0 . LT. DEL)G0TO21 

OELsD 

ALPHAs Z 

IOUTsI 

IINsI 

IF(J.LT.O)IIN«C 

LBsj 

21  CONTINUE 

IF(  IOUT . NE . 0 )G0T029 

27  CONTINUE 

0=0. 

DO  28  Is  1 (N 
LI=LT(I) 

28  Q=Q+X(LI)»(G<I)-B<LI)) 

0=  .5*0 

RETURN 

29  CONTINUE 
SIG=1 . 

IF(G(IOUT) .GT.0.)SIG=-1 . 
LIOUTsLT(IOUT) 

LIINsLIOUT 
25  CONTINUE 

SASsGC ICAC+IOUT) 
IF(K.EQ.0)G0TO31 
DO  30  1=1 ,K 

30  G(IS+I)=G(HUI)»A(IOUT.I) 

31  CONTINUE 

DO  37  I=K1 ,N 
LIsLT(I) 

IF(LI-LIOUT)32, 37,33 

32  Z=A( LI , LIOUT) 

G0T034 

33  Z=A( LIOUT, LI) 

31*  CONTINUE 

IF( K . EQ . 0 ) G0T036 
DO  35  Jsl.K 

35  Z=Z-A( I , J)#G( IS+J) 

36  G(IS*I)=Z 

37  CONTINUE 

G( IS+IOUT)=SAS 

IF(K.EQ.0)GOTO42 

G(IS+K)s-A(IOUT,K) 

IF(K.EQ.1)G0T0H2 

I=K 

00  41  11=2, K 
1=1-1 

Z=-A(IOUT,I) 

11=1*1 

DO  40  J=I1,K 
*0  Z=Z-G(ir*J)»A(J,I) 
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41  G(IS*I)»Z 

42  CONTINUE 

IF( SIG. EQ. 1 . )G0T051 
DO  50  1=1,  N 

50  G(IS+I)=-G(IS*I) 

51  CONTINUE 
IF(K.EQ.0)G0T062 
DO  61  1=1, K 

IF(G(IS+I) . EQ.O. )G0T061 
LI=LT(I) 

Jsl 


Z=BL(LI)-X(LI) 

IF(G( IS+I) . LT . 0 . ) GOT060 
J=0 


Z=BU(LI)-X(LI) 

60  CONTINUE 
Z=Z/G(IS+I) 

IF( Z.GE. ALPHA) G0T061 

ALPHA=Z 

LB= J 

IIN=I 


LIIN=LI 

61  CONTINUE 

62  CONTINUE 

X(LIOUT)=X(LIOUT)-fSlG*ALPHA 
IF(K.EQ.0)G0T071 
DO  70  1=1, K 
LI=LT( I) 

70  X(LI)=X(LI)«.ALPHA»G(IS+I) 

71  CONTINUE 


DO  72  IsKl.M 

72  G(I) =G(I)+ALPHA*G(  IAS+I) 

IF( TIN. EQ.O )G0TO90 
X(LIIN)=BL(LIIN) 

IF( LB.  EQ.O)X(LIIN)  =  BU(LIIN) 
IF( IIN.EQ. IOUT)GOTO20 
K2=(C-1 


SG=G( ID+IIN) 
I1=IIN+1 
DO  80  IsIl.N 
80  G( IV+I)=A( I , IIN) 
IF( IIN.EQ. K)G0T086 
I2=IIN+2 
S0=1  ./SG 
DO  85  IsII N,K2 
V=G(IV+I 1) 
VD=V/G(ID^I1 ) 
S1=S0-»V*VD 
RsSI/SO 

G(ID*I)=G(ID+I1)*R 
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BETAs VD/ SI 
IF(R.GT.4.)G0T0841 
DO  81  J=I2,N 

81  G(IV4j)sG(IVW)-V»A(J#I1) 
IFUl.GT.K2)G0T083 

DO  82  J*I1,K2 
J1*J+1 

82  A( J ,I)*A( J1 ,11 )«-BETA#GdV+J1) 

83  CONTINUE 
A(K,I)sBETA 
DO  84  JsKI.N 

84  A(J,I)=A(J,I1)+BETA«G(IV*J) 

G0T0849 

841  CONTINUE 
IF(I1.GT.K2)G0T0843 
DO  842  JsI1,K2 
J1sJ*1 

842  A(J,I)  =  BETA»G(IV«-J1)«.A(J1,I1)/R 

843  CONTINUE 
A(K,  I)sBETA 
DO  844  J=K1  N 

844  A(J,I)  =  BETA»GdV«-J)*A<J.I1)/R 
DO  845  J*I2,N 

845  G(IV+J)=G(IV+J)-V*A(J ,11) 

849  CONTINUE 

LT(I)sLT(I1) 

S0=S1 
11  =  12 

85  12*12*1 
SG=1 ./SI 
LT(K)=LIIN 
G( ID*K) =SG 

IF( IIN. EQ. 1 )G0T0851 
II=IIN-1 
DO  852  1*1,11 
Z=A( IIN, I) 

DO  853  J=IIN,K2 
853  A( J ,  I)  =  A(  J+1 , 1) 

852  A(K.I)*Z 
851  CONTINUE 

86  CONTINUE 

DO  87  I«K1,N 

87  G( ICAC*I)*G( ICAC+I)+SG*G( IV*I)*f2 
K1*K 

K*K2 

IIN.O 

ALPHA* 1E75 
SAS=G( ICAC+IOUT) 

IF( SAS . GT . 0 . ) ALPHA*  A  BS( G ( IOUT) > / SAS 
IF(G(IOUT) .LT.O.) G0T0898 
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J=1 

Z=X( LIOUT)-BL(LIOUT) 

G0T0899 

898  CONTINUE 
J=0 

Z : BU( L IOUT )-X(L IOUT) 

899  CONTINUE 
IF(Z.GE.ALPHA)G0T025 
ALPHA: Z 

LB:  J 

IIN:IOUT 

LIIN:LIOUT 

G0T025 

90  CONTINUE 
K2=KU1 

IF(  SIG . EQ. 1 . )G0T091 
DO  901  I:K1.N 
901  G(IAS+I)=~GUAS*I) 

91  CONTINUE 

IF(  IOUT . EQ.K1 )G0T097 
LT(IOUT)=LT(K1) 

LT(K1 )=LI0UT 
G( IAS+IOUT):G( IAS+K1 ) 

G(  ICAC+IOUT) :G(  ICAC+K 1 ) 
G(ICAC*K1):SAS 
G( IOUT)=G(K1 ) 

IF(K.EQ.0)GOTO97 
DO  92  1=1, K 
Z=A( K 1 , I) 

A(K1 , I)=A( IOUT, I) 

92  A( IOUT , I)=Z 

93  CONTINUE 
IF(K2.EQ.I0UT)G0T095 
I 1 : IOUT— 1 

DO  94  I=K2 , I 1 

94  A( IOUT ,I)=A(I,K1) 

95  CONTINUE 

IF(  IOUT. EQ. N)G0T097 

I1:IOUT*1 

DO  96  1:11 ,N 

96  A{  I , IOUT): A( I  ,K1 ) 

97  CONTINUE 
G(K1):0. 

K=K  1 

IF(  K . EQ. N)G0T027 
DO  98  I:K2,N 
Z:G(  IAS+D/SAS 
A(I,K1):Z 

98  G(  ICAC*I):G(ICAC+I)-Z*G(IAS+I) 
K1:K2 
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GOT020 

ENTRY  VE.  i 
IF(K.EQ.O> RETURN 
ID=N+N 

G(N+1 )=1  ,/G(  ID+1 ) 
IF(K.EQ.I) RETURN 
N1=K-1 

DO  111  Ial.NI 
11=1+1 

Adl.Ds-Adl.I) 

IF(I . EQ.N1 )GOT0102 
11=1+2 

DO  101  J=II,K 
Z=A(J.I) 

J1=J-1 

DO  100  L=I1,J1 

100  Z=Z+A( J,L)*A(L,I) 

101  A( J , I)=-Z 

102  CONTINUE 

AA=1 ./G( ID+I1 ) 
G(N+I1)=AA 
DO  111  Jsl.I 
Z=A( I 1 , J) *AA 

G(N+J)=G<N+J)+Z»A(I1,J) 

IFd.EQ.DGOTOIII 

J1=J+1 

DO  110  L=J1,I 

110  A(L, J)=A(L,J)+A(I1 ,L)*Z 

111  A(I1,J)=Z 
RETURN 
END 
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PROGRAM  MRRV< INPUT, OUTPUT, TAPE10) 

DIMENSION  X(16) ,EPS(16) 

COMMON/NFEST/NFEST 
COMHON/ VFO 1 E/ C( 1 50 ) 

COMMON/ VF01F/GC( 25, 50) 

C0MM0N/VF01G/T(150) 

C0MM0N/VF01 I/G2P( 325) 

COMMON/P/TF , NA , TTA(7 )  ,TA(7) ,NM,TTM(7) ,TM(7) 
COMMON/ CNSTR/ ME , MI , MC 

COMMON/TIME/ TMAX , TO , T1 , IRS , MINS , AK , ITN , IFN 
COMMON/ I NOM/INOM 
COMMON/ I PROB/IPROB 

DATA  DPR/57.29577951/ 

IPR0B=2 

IPROBsl 

IRSsI 

IRSsO 

TMAX= 1 200 . 

N=  1 1 

NA=5 

NM=5 

M=1 

MCsM 

K=1 

ME=K 

MI=M-K 

X(1)=3200./806. 48263 
DO  5  1=1,5 

TTA(  I)=r LOAT( 1-1) *0.25 
X(I+1)=20./DPR 
TTM(I)=FL0AT(I-1)*0.25 
X(I+NA+1)s60./DPR 
5  CONTINUE 
DO  1  1=1, N 
EPS(I)=1.E-4»X(I) 

1  CONTINUE 
AKMIN=1 .  E-4 
DFN=  .5 
MAXFNs 10000 
IPR1=1 
IPR2=1 
IH=2500 
M0DE= 1 

C(M*1 )*0. 136 

AK*1E60 

M2=M+M 

NN=N*(N*1 )/2 

MINSsI 

ITNsO 

IFNsO 
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PRINT  20 
20  FORMAT(IHI) 

IF(IRS.EQ.O)  GO  TO  8 

READ( 10. 16)  IRS, MINS, MODE , NFEST, ITN, IFN,(  X( I) ,1=1  ,H) ,(C(I) ,1*1  ,M2 
1  ) ,(T(I) ,1:1 ,H2) ,(G2P(I) ,1*1 ,NN) , AK.DFN 

16  F0RMAT(9X, 6(110), 135(/.9X, 5020)) 

PRINT  17.  IRS,  MINS, MODE,  NFEST,  ITN,  IFN.(XU),  1*1  ,H),(C(I),  1*1,  M2 

1  ),(T(I),I=1 ,M2) ,(G2"  I)  ,1=1, MM), AK.DFN 

17  F0RMAT(9X. 6(110), 135(/,9X,5E2«1. 1^)) 

PRINT  15 

15  FORMAT(// IX* INPUT  VALUES  OF  THETA  AND  SIGMA  USED*) 

PRINT  22, NFEST 

22  FORMAT( 1X*INPUT  VALUE  OF  HESSIAN  USED*//1X*  TOTAL  FUNCTION  • 
$*EVALUATIONS  SO  FAR* 18) 

8  CONTINUE 

CALL  SECOND(TO) 

CALL  VF01A(N,M,K,X, EPS, AKMIN.DFN.MAXFN.IPRI .IPR2.IW, MODE) 
IF(TI-TO.LT.TMAX)  GO  TO  9 
IRS:  1 

PRINT  30.TMAX 

30  F0RMAT(//1X* . TIME  LIMIT  OF*G15.8*EXCEEDED . •//) 

PRINT  31.T1-T0 

31  FORMAT( 1 X*COMPUTATION  TIME*G16.8/) 

MODEa-3 

9  REWIND  10 

WRITE(  10, 16)  IRS. MINS, MODE. NFEST,  ITN,IFN,(X(I)  ,1=1  ,F.’)  ,(C(I)  ,1*1  ,M2 
1  ),(T(I),Is1,M2) ,(G2P(I) ,1=1 ,NN) , AK.DFN 

PRINT  17.  IRS, MINS, MODE, NFEST. ITN, IFN,(X(I) ,1*1, N),(C(I), 1*1, M2 
1  ) ,(T(I) ,1=1 ,M2) ,(G2P(I) ,1=1 ,NN) , AK.DFN 

PRINT  10. NFEST 

10  F0RMAT(//5X*T0TAL  FUNCTION  EVALUATIONS* 1X16) 

END 
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SUBROUTINE  VF01B(N,M,X) 

DIMENSION  X(1) 

COMMON/NFEST/NFEST 

COMMON/ VFO 1 C/  F ,  NM ,  K ,  IS ,  MHK ,  MU 

COMMON/ VFO 1 D/G ( 50 ) 

C0MM0N/VF01E/CO50) 

COMMON/ VF0 1 F/GC<  25,50) 

COMMON/ IGAMNA/IGAMMA 
DATA  NFEST/O/ 

CALL  SG(X,F) 

IF( IGAMMA.EQ. 1 )  RETURN 

NFEST=NFEST*1 

CALL  SGX(N,X,M.NFEST) 

RETURN  $  END 


I 


SG  PAGE  1 


Si 


SUBROUTINE  SC(X.F) 

REAL  M.N.HASS 

DIMENSION  X(16).Y(9),Z(9),DY(9) 
C0MM0N/VF01E/S( 150) 

COMMON/ S$/SS( 50) 

COMMON/I NOM/ I NOM 
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COMMON/P/ TF .  NA .  TTA(  7 ) , TA(  7 ) ,  NM ,  TTM(  7 ) .  TM(  7 ) 

COMMON/ CNSTR/ME, MI, MC 


COMMON/ IGX/IGX 

COMMON/TRA J/ A , M , N . (COUNT , TIMETB e  MASS . NDE 
COMMON/ COSI/CF . VCGOR , CSI , SSI 
COMMON/IPROB/IPROB 
COMMON/ IGAMMA/ IGAMM A 

DATA  W, IGX , DPR/5 .880964 IE-2 ,0,57 .29577951/ 

IGAMMAsO 

TF=X(  1 ) 


DO  5  1=1 ,NA 
TA(  I)  =X(  I-*-1 ) 


DO  6  I-1.NM 
6  TM( I) =X( I+NA+1 ) 
TIMETB=X(  1+NA+NM+1 ) 

IF( IPR0B.EQ.2)  GO  TO  59 
TIMET9=10. 

59  CONTINUE 


NDE=6 

DELT= .004  $  NIS=250 
T=0. 

Y<1)=0.0  $  Y(2)=0.0  $  Y(3)=1. 017432502 

m!=0?T?(I)-J  Y(5)=-020943951  * 

KOUNTsO 

TIMETUTIMETB 

DO  15  1=1 , NIS 

IF(INOM.NE.I)  GO  TO  25 

IF( MOD( 1,25) .NE. 1 )  GO  TO  25 

IF(I.EQ.I)  PRINT  60 

“I;" : PSI  ■ • 

39X,* (SEC)« ,3X,*(DEG)* ,3X,#(DEG)* ,3X,*(FT)*  3X  *(FT/SEC)#  ?Y 
«  •(DEG)..3Xf.(»EO)*.3X,»(DE«»  3X  •  0K)i)  ’  ' 

CALL  DERIV(T.Y.DY) 

Z(1)=DPR«Y(1)  $  Z(2)=DPR«Y(2)  $  Z(3)=2.092642BE7*fY/?1  1  ^ 

WO.ESMT.WKO  »  Z(5)  =  OP«.T(5)ti(6).S.J(6) 
Z(7)=806.48263#Y(7)  $  Z(8)=806.48263»Y(8) 

PA=DPR»A  $  PM=DPR»M 
TIME=806 .48263»T*TF 
TEMPsVCGOR»CSI+W»CF 

CIxCF«TEMP/SQRT(TEMP**2+(VCG0R*SSI)»*2) 

PRINT  20 ,T, TIME, ( Z(J), J=1, 6), PA, PM, Z(7),Z( 8), N 
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20  FORMAT! IX , F5 . 3 »2X ,F6 . 1 .2X.F6.2 ,2X ,F6 . 3.2X ,F7 .0 ,2X,F6 .0 ,2X ,F6 .2 ,2X, 
1  F6 . 2 , 2X ,  F6 . 3  •  2X ,  F6 . 2 , 2X ,  F6 . 3 . 2X ,  F6 . 3 , 2X ,  F6 . 3) 

25  CONTINUE 

IF(KOUNT.EQ.2)  GO  TO  10 
DELT1=TIMETL/TF-T 
IF(DELT1 . LE.DELT)  GO  TO  11 

10  CALL  RUNG£(T.Y,DELT,ND£) 

IF(Y(5) .GT.-1 .)  GO  TO  50 
IGAMMAs 1 

RETURN 
50  CONTINUE 
GO  TO  15 

11  CONTINUE 

CALL  RUNGE(T.Y.DELTI.NDE) 

KOUNT=KOUNT+1 
IF C INOH.NE. 1 )  GO  TO  63 
TIMEL=806.48263*TIMETL 
IF(KOUNT.EQ.I)  PRINT  61.TIMEL 
61  FORMAT! 1H  ,*IGNITION  OCCURRED  AT  *,F10.5,*  SECONDS*) 

63  CONTINUE 

IF( DELT1 .EQ.DELT)  GO  TO  12 

DELT2=DELT-DELT1 

CALL  RUNGE( T ,  Y , DELT2 , NDE ) 

12  T IMETL=  T IMETL+O .63901697 
15  CONTINUE 

CALL  DERIV(T,Y,DY) 

TEHP= VCGOR*  CSI+W*  CF 

CI=CF*TEMP/SQRT!TEMP»*2+(VCGOR*SSI)**2> 

IF( INOM.NE. 1 )  GO  TO  26 

Z(  1 )=DPR*Y! 1 )  $  Z(2)=DPR*Y(2)  $  Z(3>=2.0926428E7*(Y!3>-1 .) 
Z(4)s25947.772*Y!4)  $  Z(5)=DPR*Y(5)  $  Z!6)=DPR*Y(6) 
Z(7)s806.48263*Y(7)  $  Z!8)s806.48263*Y(6) 

PAsDPR*A  $  PMsDPR*M 
TIME=806.48263*T*TF 

PRINT  20, T, TIME, ( Z! J) , J=1 ,6) ,PA,PM,Z(7) ,Z(8) ,N 
26  CONTINUE 
IN0M=0 

IF(  IPR0B.EQ.2)  GO  TO  72 
Fs-Y(2) 

SS(1)s!Y!3>-1. 0047786464)70. 0047786464 
SS(2)=Y(7) 

S S(3)=Y(8) 

GO  TO  74 
72  CONTINUE 
FsCI 

SS(1)=(Y(3)-1. 029054170)/. 029054170 
SS(2)=Y(4)/. 90566542-1 . 

SS(3)=Y(5) 

SS(4)sY(7) 

SS(5)*Y(8) 
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74  CONTINUE 

IF< IGX. EQ. 1 )  GO  TO  40 
DO  35  1=1 ,MC 
S(I)=SS(I) 

35  CONTINUE 


PRINT  28,Ff(S( J) ,Jsl ,MC) 

28  F0RMAT(1X./,1X.6F20.15./) 
40  CONTINUE 


RETURN 

END 
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SUBROUTINE  RUNGE  (T.X.DELT.N) 

DIMENSION  X(10) ,DX< 10) ,DELX( 10,3) .XV( 10) 

CALL  DERIV  (T,X,DX) 

T2  s  T  ♦  DELT/2.0 
DO  100  I  s  1,N 
DELX(I.I)  =  DX( I)*DELT 
100  XV(I)  *  X( I)  ♦  DELX<I,1)/2.0 
CALL  DERIV  (T2.XV.DX) 

DO  200  Isl  , N 
DELX<I,2)  =  DX(I)  *DELT 
200  XV(I)  =  X(I)  ♦  DELX( I ,2)/2.0 
CALL  DERIV  (T2.XV.DX) 

DO  300  1=1 ,N 
DELX(If 3)  =  DX(I)»DELT 
300  XV(I)  =  XU)  +  DELX( 1,3) 

T  =  T  DELT 

CALL  DERIV  (T.XV.DX) 

DO  400  1=1 ,N 

400  X( I)  =  X(I)  ♦  (DELX(I.I)  ♦  DX(I)»DELT  ♦  2 .0*(DELX(I,2)  ♦  DELXU.3) 
1))/6.0 
RETURN 
END 
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SUBROUTINE  DERIV(TT,Y,DY) 

REAL  L.M.N.NLIH.MASS 
DIMENSION  Y(9>,DY<9> 

COMMON/P/TTF , NA, TTA( 7 )  ,TA(7> ,HH,TTM(7>  ,TH(7> 

COMMON/ TR A J/ A , M , N , KOUNT , TIMETB , MASS , NDE 
COMMON/ COSI/CF , VCGOR , CS , SS 

DATA  TW.W2, ALIM, NLIH/. 11761928. 3. A585739E-9.. 69813170, A. 5/ 

TsY( 1 )  $  FsY<2)  $  R*Y(3>  8  V«Y<4)  8  GsY(5>  8  SsY(6) 

SF=SIN(F)  $  CF=COS(F> 

SGtSIN(Gl  $  CG=COS(G) 

SStSIN(S)  $  CSsCOS(S) 

TF=SF/CF  8  TGsSG/CG 
CALL  SLIN1(TTtA,TTA,TA,NA) 

CALL  SLINKTT.M.TTM.TM.NH) 

SAsSINCA)  8  CA=COS( A)  8  SMrSIN(M)  8  CMsCOS(M) 

VCG=V*CG 

VCGOR=VCC/R 

CALL  AERO^S,V,A,SA.CA,D,L> 

GRs 1 ./R**2 

IF(KOUNT.GT.O)  GO  TO  5 
MASS: 1 . 

THF=0. 

GO  TO  15 

5  IF(KOUNT.GT.I)  GO  TO  10 

MASS:  1 . - . 83533982* ( TT*TTF-TIMETB) 

THR=. 30555556 
GO  TO  15 

10  MASS:  .*<6620367 
THR=0. 

15  CONTINUE 
N=L/MASS 

TMDOM= ( THR* CA-D) /MASS 
TPLOM=(THR*SA-*-L)/MASS 
TDOT  s  V  CGOR*  CS/ CF 
FDOTsVCGOR*SS 
RDOT=V*SG 

VDOTsTMDOM-GR*SG>W2*R*CF* ( SG*CF-CG*SF*SS) 

CDOTs  TPLOM*CM/  V-G  R*  CG/V+VCGOR+TW*  CF*  CS*W2*  (  R*  CF/ V )  •  (  CG*  CF+SG»SF*SS 

1  ) 

SDOT  s  T  PLOM*SM/  V  CG-VCGOR*  TF*  CS+TW*  ( TG*SS*  CF-SF)  -  ( W2/V  CGOR)  *SF*  CF*  CS 

DY( 1 )=TDOT*TTF  8  DY(2)=FDOT*TTF  8  DY(3)=RDOT*TTF 

DY( A )=VDOT*TTF  8  DY(5)=GDOT*TTF  8  DY(6)sSDOT*TTF 

IF(NDE. EQ.6)  GO  TO  20 

P: 1 .-A/ ALIM 

IFCP.GT.O.)  PsO. 

DY(7)=-P*P 

PsN 

IF(P.GT.O.)  P=0. 

DY(8)=-P*P 

DY(7)=DY(7)#TTF  $  DY(8)=DY(8)*TTF 
20  CONTINUE 

RETURN  8  END 
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SUBROUT I NE  AEROC  R , V , AL , SAL , CAL . D . L) 

REAL  L.M 

DIMENSION  TM(5) ,TAA(5) ,TBB<5> ,TCC(5) ,TA1(5) ,TA2(5> .TB1<5) ,TB2<5> ,T 
ACK5)  ,TC2(5>  ,TD1<5)  ,TD2<5)  ,TD3<5) 

COMMON/IPROB/ I PROS 
DATA  TM/. 2, 1.2, 5.. 10., 20./ 

DATA  TAA/0.,0.,0.,2.4E-12,1.38E-12/ 

DATA  TBB/7 .8E-8 ,7 .7E-8 .5.6E-8 .-7 . 1 IE-7 .-3 .81E-7/ 

DATA  TCC/ .01 IN , .0076 , .0028 , .0578 , .0297/ 

DATA  TA1/-1 .OE-4 ,-7.8lE-5,1 .09E-5.1 .41E-5.1 .33E-5/ 

DATA  TB1/-2 . 19E-3 ,-1 .25E-3 .-9 .62E-4 ,-9 .0E-4 ,-8 . 19E-4/ 

DATA  TA2/- 1 . OE-4 , -7 . 8 1 E-5 , 4 . 86E-6 , -6 . 6E-6 , -6 . 25  E-6/ 

DATA  TB2/-2 . 19E-3 ,-1 . 25E-3 ,-1.1 3E-4 , 3 . 49E-4 , 3 .58E-4/ 

DATA  TCI/ .035 . .076 , .0324 , .0261 , .0243/ 

DATA  TC2/ .035 . .076 , .0204 , .011 4 , .0 1 05/ 

DATA  TD1 / 1 . 33E-4 ,-7 . 8 1 E-5 , 2 . 94E-4 , 4 . 38E-4 , 4 . 38E-4/ 

DATA  TD2/2 . 68E-2 , 2 . 90E-2 , 1 . 4 1 E-2 , 6 . 50E-3 , 6 . 50E-3/ 

DATA  TD3/- . 030 , - . 030 , - . 035 , - . 035 . - . 035/ 

DATA  DPR.S/57. 29577951. 125. 84/ 

Hs20926428.«(R-1.) 

VEL=25947 .772*V 
ALFA=DPR* AL 

CALL  ATM62(H.TAU. SIGMA) 

M= VEL/( 65 . 770*SQRT< TAU) ) 

CALL  SLINK M.AA,TM,TAA,5) 

CALL  SLIN1(M,BB,TM,TBB,5) 

CALL  SLIN 1 ( M , CC , TM , TCC , 5 ) 

IF( ALFA.GT. 16,)G0  TO  10 
CALL  SLINKM.A.TM.TAI  ,5) 

CALL  SLIN1(M,B,TM.TB1.5) 

CALL  SLINK M.C.TM.TC1 ,5) 

GO  TO  20 

10  CALL  SLINKM,A,TM,TA2,5> 

CALL  SLINKM.B.TM.TB2.5) 

CALL  SLINKM,C,TM.TC2,5> 

20  CONTINUE 

CALL  SLINKM.D1  ,TM,TD1 ,5) 

CALL  SLINK M,D2,TMfTD2,5) 

CALL  SLINKM.D3.TM.TD3. 5) 

CASF=AA*H*H+BB*Hi-CC 
CAPRs  A*  ALFA* ALFA+B* ALFA+C 
CArCASF+CAPR 

CN= D 1 • ALFA*  ALFA+D2*ALFA+D3 
CL«CN*CAL-CA*SAL 
CDsCA*CAUCN*SAL 
IF(IPR0B.EQ.2>  GO  TO  30 
BB* 1 5268 .635*SIGMA*V**2 
GO  TO  40 

30  BB=932.34367*SIGMA*V**2 

40  CONTINUE 
D*BB*CD 
LsBB*CL 
RETURN 
END 
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SUBROUTINE  ATM62(H,T,S) 

REAL  LHC 8) 

DIMENSION  HH(8) ,TMC8) , CSC 8 ) 

DATA  HH/O . , 36089 .  ,65617 . .  10*1987 .  ,15*1199.  •  170604 .  ,2001 31 .  ,259186 ./ 
DATA  TH/288.15,216.65.216.65,228.65,270.65,270.65.252.65,180.65/ 
DATA  LM/-1.9812E-03,0..+3.048E-04,8.5344E-04,0.,-6.096E-04, 

1  -1.2192E-O3.0./ 

DATA  CS/3.401 824655257E-1 1 , 1 .683376997 1 49E+00 , 
19.817858914969E+80,1.506414967722E>29,4.394749884481E-O1, 
24.717851 690435E-43 , 1 . 56279374065 1 E-22 ,5 .028946984 1 09E+01 / 

DATA  GOR/ .010413309/ 

IFCH.GE.O.)  GO  TO  5 
Ts288. 15-1 .9812E-3*H 
S=1 .-2.926291 3E-5*H 
RETURN 

5  IFCH.LE. 750000.)  CO  TO  10 
T=180.65 
S=8 . 1 1 1816E-18 
RETURN 
10  CONTINUE 
DO  20  Is1,8 
20  IF(H.GE.HM( I) )M=I 

T=TM( M) ♦LHC  M) • ( H-HMC M) ) 

IF(LMCH) .EQ.O. >30.40 
30  SsCSC M) *EXPC-GOR*H/T) 

RETURN 

40  S=CS( M) *T#*C-1 .-GOR/LMC M) ) 

RETURN 

END 
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SUBROUTINE  SLIHKX,  Y.  AX,  AY.  N> 

C  SINGLE  LINEAR  INTERPOLATION  SUBROUTINE 
DIHENS ION  AX(N),  AY(N) 

IF(N.EQ.I)  GO  TO  30 
IF(X.LE.AXO))  00  TO  20 
IF(X.GT.AX(N))GO  TO  30 
DO  10  2«1,N 
K  «  I 

Z  «  X  -  AX(I) 

IF( 2.LE.O.)GO  TO  11 

10  CONTINUE 

11  J«K-1 

S  s  ( AY(K>  -  AY( J) >/( AX(K)  -  AX(J)> 

Y  •  AY( J)  ♦  S«(  X  -  AX( J) ) 

GO  TO  100 
20  YsAY(1) 

GO  TO  100 
30  Y*AY(N) 

100  RETURN 
END 
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SUBROUTINE  SGX(N ,X,H,NFEST) 

DIMENSION  X( 16)  ,XF( 16) ,XB( 16) ,CF( 16) ,CB( 16) 
COMMON/VFOI D/G(50) 

COMMON/VF01F/GC(25,50) 

COHMON/SS/SSC50) 

COMMON/ IGX/IGX 
DATA  EPSP/1.E-V 
IGXrl 

DO  15  1*1, N 
XF(I)sX(I) 

15  XB(I)=X(I) 

DO  20  1*1. N 
DX:ABS( EPSP*X( I) ) 

IF( ABS(X(I) ) . LE.EPSP)  DX=EPSP“2 
XF( I)=X(I)+DX 
CALL  SG(XF.FF) 

DO  16  Jsl.M 

16  CF( J)sSS( J) 

XF(X)«X(X) 

XB( I)=X( I)-DX 
CALL  SG(XB.FB) 

NFEST=NFEST*2 
DO  17  J=1,M 

17  CB( J)*SS( J) 

XB(I)sX(I) 

G(I)=0.5»(FF-FB)/DX 
DO  20  J=1.M 

GC( I,J)=.5*(CF(J)-CB(J)) /DX 
20  CONTINUE 
IGX=0 
RETURN 
END 


APPENDIX  D 

GLOSSARY  OF  VARIABLES  IN  USER  CODE 


102 


PR  JGRAM  MRRV 


Jyiu\ 

nef in  It  ion 

AK 

Real 

Measure  of  constraint  convergence  rate 

AKMIN 

Real 

Defines  convergence  for  constraints 

C(I) 

Real 

Constraint  residuals  and  constraint  scale  factors 

DFN 

Real 

Expected  decrease  in  performance  index 

DPR 

Real 

Degrees  per  radian 

EPS(I) 

Real 

Defines  convergence  in  VA09A 

GC(I) 

Real 

Derivatives  of  the  constraints 

G2P(I) 

Real 

Second  derivative  matrix  In  factored  form 

I 

Integer 

Counter 

IFN 

Integer 

Number  of  function  evaluations  in  VA09A 

INOM 

Integer 

Trajectory  print  flag.  Print  occurs  if  INOM  •  1 

IPROB 

Integer 

Problem  flag.  »  1,  reentry;  ■  2,  plane  change 

IPR1 

Integer 

Print  flag  in  VF01A 

IPR2 

Integer 

Print  flag  in  VA09A 

IRS 

Integer 

Restart  flag.  *  0,  first  run;  -  1,  restart 

ITN 

Integer 

Number  of  iterations  in  VA09A 

IW 

Integer 

Number  of  elements  in  work  space  W(I) 

K 

Integer 

Number  of  equality  constraints 

M 

Integer 

Number  of  constraints 

MAXFN 

Integer 

Maximum  number  of  function  evaluations  in  VA09A 

MC 

Integer 

-  M 

ME 

Integer 

Number  of  equality  constraints 

Ml 

Integer 

Number  of  inequality  constraints 

MINS 


Integer 


Number  oE  iterations  in  VF01A 


PROGRAM  MRRV,  Cont'd 


Variable 

IZ£i 

Definition 

MODE 

Integer 

Flag  which  indicates  how  6,  o,  and  G2P  are  input 

M2 

Integer 

-  2M 

N 

Integer 

Number  of  parameters 

NA 

Integer 

Number  of  parameters  in  angle  of  attack  table 

NFEST 

Integer 

Number  of  function  evaluations 

NM 

Integer 

Number  of  parameters  in  bank  angle  table 

NN 

Integer 

Number  of  elements  in  G2F 

T(I) 

Real 

Array  where  &  and  a  are  stored 

TA(I) 

Real 

Node  ordinates  for  angle  of  attack  table 

TF 

Real 

Pinal  time 

TM(I) 

Real 

Node  ordinates  for  roll  angle  table 

TMAX 

Real 

Internal  time  limit 

TTA(I) 

Real 

Node  abscissas  for  angle  of  attack  table 

TTM(l) 

Real 

Node  abscissas  for  roll  angle  table 

TO 

Real 

Initial  value  of  internal  time 

T1 

Real 

Current  value  of  internal  time 

X(I) 

Real 

Initial  values  of  unknown  parameters 

SUBROUTINE 

VF01B 

Variable 

Im. 

Definition 

C(I) 

Real 

Constraints  and  constraint  scale  factors 

F 

Real 

Performance  index 

C(I) 

Real 

Derivatives  of  F  with  respect  to  X(I) 

GC  (I , J) 

Real 

Derivatives  of  C(J)  with  respect  to  X(I) 

IGAMMA 

Integer 

Flag;  •  0,  y  >  -60  deg  ;  ■  1,  y  <  -60  deg 
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SUBROUTINE  VFOIB.  Cont'd. 


Variable 

Type 

Definition 

IS 

Integer 

Used  in  optimization  code 

K 

Integer 

Number  of  equality  constraints 

M 

Integer 

Number  of  constraints 

MM 

Integer 

Used  in  optimization  code 

MMK 

Integer 

Used  in  optimization  code 

N 

Integer 

Number  of  parameters 

NFEST 

Integer 

Number  of  function  evaluations 

NU 

Integer 

Used  in  optimization  code 

X(I) 

Real 

Unknown  parameters 

SUBROUTINE  SG 

Variable 

Type 

Definition 

A 

Real 

Angle  of  attack  (rad) 

CF 

Real 

Cos 

Cl 

Real 

Cos  i 

CSI 

Real 

Cos  $ 

DELT 

Real 

Integration  step  size 

DELT1 

Real 

Partial  integration  step  before  engine  is  started 
or  shut  off 

DELT2 

Real 

Partial  Integration  step  after  engine  is  started 
or  shut  off 

DPR 

Real 

Degrees  per  radian 

DY(I) 

Real 

Right-hand  sides  of  the  differential  equations 
of  motion 

F 

Real 

Performance  index 
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SUBROUTINE  SG,  Cont’d. 

Variable 

Type 

Definition 

I 

Integer 

Counter 

IGAMMA 

Integer 

Flag;  -  0,  Y  >  -60  deg  ;  »  1,  y  _>  -60  deg 

IGX 

Integer 

Flag;  ■  1,  prevents  update  of  constraint  residuals 
when  derivatives  are  being  computed 

INOM 

Integer 

Trajectory  print  flag;  •  1,  causes  trajectory  to  be 
printed 

IPROB 

Integer 

Problem  flag;  ■  1,  reentry;  ■  2,  plane  change 

J 

Integer 

Counter 

KOUNT 

Integer 

Engine  flag;  ■  0,  engine  off;  ■  1,  engine  on;  ■  2, 
engine  off  again 

M 

Real 

Bank  angle  (rad) 

MASS 

Real 

Vehicle  mass  (nondlmenslonal) 

MC 

Integer 

Number  of  constraints 

ME 

Integer 

Number  of  equality  constraints 

MI 

Integer 

Number  of  inequality  constraints 

N 

Real 

Load  factor 

NA 

Integer 

Number  of  parameters  in  angle  of  attack  table 

NDE 

Integer 

Number  of  differential  equations 

NIS 

Integer 

Number  of  integration  steps 

NM 

Integer 

Number  of  parameters  in  bank  angle  table 

PA 

Real 

Angle  of  attack  (deg) 

PM 

Real 

Bank  angle  (deg) 

S(I) 

Real 

Constraints  (=  C(I)) 

SS(I) 

Real 

Temporary  values  of  constraints 

SSI 

Real 

Sin  i p 

T 

Real 

Time  (normalized) 
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SUBROUTINE  SG. 

Cont 'd. 

Variable 

Type 

Definition 

TA(I) 

Real 

Node  ordinates  in  angle  of  attack  table  (rad) 

TEMP 

Real 

Temporary  computation 

TP 

Real 

Final  time  (nondimenslonal) 

TIME 

Real 

Time  (sec) 

T1MEL 

Real 

Ignition  time  or  burnout  time  (sec) 

TIMETB 

Real 

Ignition  time  (nondimenslonal) 

TIMETL 

Real 

Ignition  time  or  burnout  time  (nondimenslonal) 

TM(I) 

Real 

Node  ordinates  of  bank  angle  table  (rad) 

TTA(I) 

Real 

Node  abscissas  of  angle  of  attack  table  (nondimenslonal) 

TTM(I) 

Real 

Node  abscissas  of  bank  angle  table  (nondimenslonal) 

VCGOR 

Real 

V  cosf^r  (nondimenslonal) 

W 

Real 

Angular  velocity  of  earth  (nondimenslonal) 

X(I) 

Real 

Unknown  parameters 

Y(I) 

Real 

State  variables  (nondimenslonal) 

Z(I) 

Real 

State  variables  (dimensional) 

SUBROUTINE  RUNGE 

Variable 

Type 

Definition 

DELT 

Real 

Stepsize 

DELX(I,J) 

Real 

Increment  in  X(I) 

DX(I) 

Real 

Derivative  of  X(I) 

l 

Integer 

Counter 

N 

Integer 

Number  of  differential  equations 
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SUBROUTINE  RUNCE,  Cont'd. 


Variable 

Type 

Definition 

T 

Real 

Independent  variable 

T2 

Real 

Intermediate  value  of  independent  variable 

X(I) 

Real 

Dependent  variables 

XV(I) 

Real 

Intermediate  value  of  X(I) 

SUBROUTINE 

DERIV 

Variable 

Type 

Definition 

A 

Real 

Angle  of  attack  (rad) 

ALIM 

Real 

Angle  of  attack  limit  (rad) 

CA 

Real 

Cos  o 

CF 

Real 

Cos  $ 

CG 

Real 

Cos  y 

CM 

Real 

Cos  p 

CS 

Real 

Cos  V 

D 

Real 

Drag  (nondimensional) 

DY(I) 

Real 

Right-hand  sides  of  the  equations  of  motion 

F 

Real 

Latitude,  $  (rad) 

FOOT 

Real 

• 

<t>  (nondimensional) 

GR 

Real 

Acceleration  of  gravity  (nondimensional) 

KOUNT 

Real 

Engine  flag;  «  0,  off;  ■  1,  on;  ■  2,  off  again 

L 

Real 

Lift  (nondimensional) 

M 

Real 

Bank  angle  (rad) 

i 

s 
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SUBROUTINE  DERIV,  Cont'd. 


Vi l  r  1  uli  1 1' 

Definition 

HASS 

Real 

Maas  (nondimenalonal) 

N 

Real 

Load  factor 

NA 

Integer 

Number  of  nodes  In  angle  of  attack  table 

NDE 

Integer 

Number  of  differential  equations 

NLIM 

Real 

Load  factor  limit 

NM 

Integer 

Number  of  nodes  in  bank  angle  table 

P 

Real 

Temporary  computation 

R 

Real 

Distance  from  center  of  earth  (nond linens ional) 

ROOT 

Real 

r  (nondimensional) 

S 

Real 

Heading  angle,  \J>  (rad) 

SA 

Real 

Sin  a 

SDOT 

Real 

• 

<j/  (nondimen s ional) 

SF 

Real 

Sin  4> 

SG 

Real 

Sin  y 

SM 

Real 

Sin  u 

SS 

Real 

Sin  iji 

T 

Real 

Longitude,  9  (rad) 

TA(I) 

Real 

Node  ordinates  in  angle  of  attack  table  (rad) 

TDOT 

Real 

9  (nondimensional) 

TF 

Real 

Tan  $ 

TG 

Real 

Tan  y 

THR 

Real 

Thrust  (nondimensional) 

TIMETB 

Real 

Ignition  time  (nondimensional) 

TM(I) 

Real 

Node  ordinates  in  bank  angle  table  (rad) 
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SUBROUTINE  DERIV,  Cont'd. 


Variable 

Definition 

TIM.OM 

Kf.i  1 

(T  air  (i  +  L)/m 

TMDOM 

Real 

(T  coa  a  -  D)/m 

TT 

Real 

Normalized  time,  t 

TTA(I) 

Real 

Node  abscissas  in  angle  of  attack  table 

TTF 

Real 

Final  time  (nondimensional) 

TTM(I) 

Real 

Node  abscissas  in  bank  angle  table 

TO 

Real 

2<i)  (nondimensional) 

V 

Real 

Velocity  (nondimensional) 

VCG 

Real 

V  cos  y 

vcgor 

Real 

V  cos  y/r 

VDOT 

Real 

V  (nondimensional) 

W2 

Real 

w2  (nondimensional) 

Y(I) 

Real 

Nondimensional  states 

SUBROUTINE 

_AER0 

Variable 

Jype 

Definition 

AL 

Real 

Angle  of  attack  (rad) 

ALFA 

Real 

Angle  of  attack  (deg) 

A1 

Real 

Coefficient  in  C. 

SF 

A2 

Real 

Coefficient  in  C. 

Kr, 

A3 

Real 

r  K 

Coefficient  in  C„ 

B1 

Real 

Coefficient  in  C 

A_  _ 

SF 

P2 

Real 

Coefficient  in  C. 

\>R 

B3 

Real 

rK 

Coefficient  in  C„ 

CA 

Real 

Axial  force  coefficient 
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SUBROUTINE  AERO. 

Cent 'd. 

Vnf  1  <it<  1  »• 

Oof f nit  ton 

CAL 

Real 

Cos  a 

CAPR 

Real 

Axial  force  coefficient,  pressure 

CASF 

Real 

Axial  force  coefficient,  skin  friction 

CD 

Real 

Drag  coefficient 

a 

Real 

Lift  coefficient 

CN 

Real 

Normal  force  coefficient 

Cl 

Real 

Coefficient  in  C, 

SF 

C2 

Real 

^efficient  in  CAp^ 

C3 

Real 

Coefficient  in  C^ 

D 

Real 

Drag  (nondimensional) 

DPR 

Real 

Degrees  per  radian 

H 

Real 

Altitude  (ft) 

IPROB 

Integer 

Problem  flag;  -  1,  reentry;  -  2,  plane  change 

L 

Real 

Lift  (nondimenaional) 

M 

Real 

Mach  number 

R 

Real 

Distance  from  center  of  earth  (nondimenaional) 

SAL 

Real 

Sin  a 

SIGMA 

Real 

Density  ratio 

TAU 

Real 

Absolute  temperature  (deg  K) 

TA1(I) 

Real 

Table  for  A1 

TA3(I) 

Real 

Table  for  A3 

TBl(I) 

Real 

Table  for  Si 

TB3(I) 

Real 

Table  for  S3 

TCI (I) 

Real 

Table  for  Cl 

Ill 


SUBROUTINE  AERO ,  font 'd. 


Variable 

TyP.e 

TC3 

Real 

TM(I) 

Real 

T1A2(I) 

Real 

TxB2(I) 

Real 

T1C2(I) 

Real 

T2A2(I) 

Real 

T2B2(I) 

Real 

T2C2(I) 

Real 

V2L 

Real 

u 

Real 

Definition 

Table  for  C3 
Table  for  Mach  number 
Table  for  A2,  a  <_  16  deg 
Table  for  B2,  o  _<  16  deg 
Table  for  C2,  at  <_  16  deg 
Table  for  A2,  a  >  16  deg 
Table  for  B2,  a  >  16  deg 
Table  for  C2,  a  >  16  deg 
Velocity  (ft/sec) 

Velocity  (nondimenslonal) 


SUBROUTINE  ATM62 


Variable 

TyPe 

Definition 

CS(I) 

Real 

Coefficients  in  density  formula  (nondimenslonal) 

GOR 

Real 

*s/R 

H 

Real 

Altitude  (ft) 

HM(I) 

Real 

Altitude  at  beginning  of  layer  m  (ft) 

I 

Integer 

Counter  1 

i 

Temperature  gradient  in  layer  m  (deg  K/ft) 

LM 

Real 

M 

Real 

Denotes  layer 

S 

Real 

Density  ratio 

T 

Real 

Absolute  temperature  (deg  K) 

TM 

Real 

Temperature  at  beginning  of  each  layer  (deg  K)  J 
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SUBROUTINE 

SLIN1 

Variable 

Type 

Definition 

AX(I) 

Real 

Node  abscissae 

AY  (I) 

Real 

Node  ordinates 

I 

Integer 

Counter 

J 

Integer 

Counter 

K 

Integer 

Counter 

N 

Integer 

Number  of  points  in  table 

S 

Real 

Slope 

X 

Real 

Abscissa 

Y 

Real 

Ordinate 

2 

Real 

Temporary  computation 
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